ここ、数ヶ月、本連載がストップしてしまいました。大変、失礼をいたしました。
体調も、ほぼ、元通りになり、また、数式処理ソフト DERIVE（デライブ）を続けていきたいと思います。よろしく、お願いいたします。
なお、33項では、34項で「2階の常微分方程式のユーティリティファイル」を取り上げると書きましたが、その前に連立微分方程式を取り上げました。
ご了承下さい。
惑星の運動でもそうですが、微分方程式が1つではなく、複数個の変数に対する連立微分方程式の形をとることが、しばしば、あります。
また、1つの変数に対する、2階の微分方程式を便宜上、2つの変数に対する連立1階微分方程式として、解くことの方が有効な場合もあります。たとえば、変数として位置とその時間微分である速度との組み合わせを考えてみてください。
この項では、そのような例を取り上げて参りましょう。
再び、線形バネで束縛された質量mの質点の直線上の運動を考えます。
質点の位置（の変位）をx（t）で、速度をv（t）で表しましょう。33項のところで見たように、運動方程式は、おなじみの、mx''＝－kxです。
2階微分を''で表しています。初期条件は、t＝0で、x＝0、x’＝v0としましょう。v＝x'であるので、この2階の微分方程式を次のように1階の連立方程式に直して考えます。（33項のところと異なり、速度の初期値をv0で表しています）
［mv’＝－kx， x’＝v］
上記を演算子法で解きましょう。xの像関数（ラプラス変換）をF（s）、vの像関数をG（s）とします。
すると、与式は、初期条件を考慮して、［m（sG－v0）＝－kF，sF＝G］となります。DERIVEでこの連立方程式を解くと、
[F = mv0/(k + ms^2)，G = msv0/(k + ms^2)]となりました。ラプラス逆変換のためのユーザ定義ファイルを開いておけば、これより、
INVLaplace（式、s、t）として、直ちに、x＝√mv0SIN(√kt/√m)/√k、v＝v0COS(√kt/√m) が得られます。ここで、式とは、FまたはGのことです。
重力の働かない空間でのロケットの運動（直線上を運動します）を取り上げてみます。
燃料の質量をm（t）、初期値をm0、ロケットの燃料以外の部分の質量をMとします。
また、燃料の噴出速度（ロケットに対する相対的な速さ）をu、単位時間に噴出する燃料の質量をkとして、uとkは、簡単なために一定としましょう。
さらに、ロケットの速度をv（t）で表します。この初期値は、ゼロとしておきましょう。
このとき、ロケットには、他から力が働かないため、ロケットと噴出した分の燃料との運動量の変化分は、ゼロとなります。
すると、静止した固定座標系から見たとき、時刻tのときのロケットの運動量は、（M＋m）v、（進行方向を正とします）。
t＋dtのとき、ロケットと噴出した分の燃料との運動量の計は、（M＋m－kdt）（v＋dv）－k（u－v）dt、
その差をゼロと置いて、（M＋m）dv－kvdt－kudt＋kvdt＝0、
よって、解くべき方程式は、dv／dt＝ku／（M＋m）となります。ここで、mは、m＝m0－ktです。T＝m0／kとなるt＝Tで燃料はなくなります。
これより、v（t）＝uLN((M + m0)/(M + m0 － kt))
最終到達速度の大きさは、uLN（1＋m0／M）です。これは、uに比例して、kには依存しません。最終到達速度は、uに正比例するため、噴射速度の大きさが重要なポイントとなります。
なお、kが大きいと、最終到達速度に達するまでの時間Tが短縮できます。
地表からと言っても、他の惑星でも同様です。ただし、簡単なために地球の大気による空気抵抗及び自転の影響を無視します。
このとき、（M＋m）dv／dt＝ku－G（M＋m）Me／（R＋z）^2
ただし、垂直方向にz軸を取り（上方を正）ます。ここで、G：万有引力定数、Me：地球の質量、R：地球の半径（地球は真球と考えます）です。
また、m＝m0－kt、v＝dz／dtの関係があります。
zがRに比して小さいときを考えてみましょう。このとき、dv／dt＝ku／（M＋m0－kt）－g となります。ここで、g＝GMe／R^2で地表における重力加速度です。
これより、v＝－ uLN(－M + kt － m0) － gt＋C、Cは積分定数です。初速度を0とすれば、v＝ uLN((M + m0)/(M
+ m0－kt)) － gt 。
ただし、この式は、t＝0～m0/kまでで、それ以降は、v＝uLN((M + m0)/M) － gm0/k－g（t－m0/k） （自由落下）となるわけです。
少なくともロケットが上方に打ち上がるためには、t＝0で正のvが必要となります。この条件から、k/(M+m0)>=g/uが得られ、整理すると、uk>=（M+m0）gがその必要条件となります。このとき、t＝m0/kでも正の速度を維持していますので、それから、しばらくは、上方に運動を続けて、最高高度に達し、その後に自由落下することになります。最高高度に達する時刻は、t
= uLN((M + m0)/M)/g です。
t＝0～m0/kの間で、vの式をtで積分すると、z＝－ u(M － kt + m0)LN(1/(M － kt + m0))/k － u(M
－ kt0 + m0)LN(M + m0)/k － gt0^2/2 + tu が得られます。
これより、t＝m0/kの時の高さは、－ MuLN((M + m0)/M)/k － gm0^2/(2k^2) + m0u/k です。
上の結果を使って、t>m0/kの時のvの式を積分すると、z＝－ u(M － kt + m0)LN(M + m0)/k + u(M － kt
+ m0)LN(M)/k － gt^2/2 + m0u/k 。
そして、到達できる、最高高度は、上式に最高高度に達する時刻を代入して、
zmax＝u^2LN(M + m0)^2/(2g) － LN(M + m0)(u^2LN(M)/g + u(M + m0)/k) + u^2LN(M)^2/(2g)
+ u(M + m0)LN(M)/k + m0u/k と求められます。
最終更新日 2008/3/11