数式処理ソフト　DERIVE（デライブ） de ドライブ

３２．1階常微分方程式の解法

１．前回の補足

「数式処理ソフト DERIVE（デライブ）の３１．運動方程式の非線形項がある場合の解は、どうなっているのかしら？」

「おう、そうじゃのう。バネの力が線形の時は、三角関数が現われたが、3次の非線形項がある場合は、楕円関数が現われるのじゃよ」

「楕円関数って、何？」

「元々は、楕円の周の長さ（弧長）を計算しようというところから研究が始まったので、こういう名前が付いているが、三角関数が三角形と強い関係があるのに比べて、楕円関数と楕円との間にそれほど深い関係はないとのことじゃ。
いずれ、別の場所で取り上げようかの。
当面、三角関数を複雑にしたような周期関数が楕円関数で、それを使えば、解を明示的に表せるということを覚えておいて欲しい」

「前回は、速度と位置とのグラフは、楕円だったけど、これは、どうなるのかしら？」

「速度（正確には、運動量）と位置座標のグラフを一般に相空間（平面）というのじゃが、今回の例で、相平面上の軌道（楕円に似ている）が囲む図形の面積を計算しようとすると、√（4次整式）の積分が現われる。
√（3次整式）、1／√（3次整式）、1／√（4次整式）の積分は、一般には、楕円積分で表されるのじゃ。
（その具体的な形は、たとえば、「数学辞典　第3版」（岩波書店）の付録に表が載っているので参照）
なお、以前、取り上げたように、DERIVEにも楕円積分がユーティリティファイルに定義されているので、楕円積分は、直ちに計算可能じゃが、これについても楕円関数のところで扱おうかの」

２．DERIVEでの解法

「さて、本題じゃ。
前回のように、数学的に式を変形して微分方程式を解くことも、もちろん、よいが、よく知られたタイプの微分方程式は、簡単に解くためのユーティリティファイルがDERIVEには、あるのじゃな。
ここでは、1階の微分方程式の一部を紹介しよう。（主としてDERIVEのヘルプと「常微分方程式）の解法」（培風館）」を参考にした。）
まずは、ファイルから「開く」で「FirstOrderODEs.mth」を開く。yを変数xの関数とする（以下、同じ）」

(0) 基本的な解法 DSOLVE1(p, q, x, y, x0, y0) または DSOLVE1_GEN(p, q, x, y, c)

p(x, y) + q(x, y)y' = 0なる微分方程式の解法。x＝x0のときにy＝y0。
なお、類似のものに、DSOLVE1_GENがあるが、これは、一般解として定数cが必要な場合に使う。（以下、_GENが付いているものは同様）
さて、たとえば、y^2ｙ’/x－y^3=0では、DSOLVE1(－xy^3, y^2, x, y, x0, y0, a_)として、
y = y0e^(x^2/2 － x0^2/2)を得る。

（1）y' + p(x)y = q(x)の場合－－－LINEAR1(p, q, x, y, x0, y0) または LINEAR1_GEN(p, q, x, y, c)

たとえば、p(x)＝－x、q(x)＝2x、y(0)＝2のときは、LINEAR1(－x, 2x, x, y, 0, 2)と数式入力行に入力して、計算させるとy = 4e^(x^2/2) － 2が得られる。
上記の例では、LINEAR1_GEN(－x, 2*x, x, y, c)と入力して計算させると、y = ce^(x^2/2) － 2が得られる。

（2）y'＝p(x)q(y)の場合－－－－－SEPARABLE(p, q, x, y, x0, y0) または SEPARABLE_GEN(p, q, x, y, c)

変数分離形の場合である。
たとえば、p(x)＝x、q(y)＝y^2では、SEPARABLE(x, y^2, x, y, x0, y0)と入力する。計算してyについて解けば、
y = － 2y0/(y0x^2 － x0^2y0 － 2)となる。ここで、x0、y0は、初期値である。

（3）y'＝r（x，y）の場合－－－－－HOMOGENEOUS(r, x, y, x0, y0) または HOMOGENEOUS_GEN(r, x, y, c)

HOMOGENEOUSとは、日本語では、「同次」と訳されている。y/xを変数とするように変形できるもの。
たとえば、（x+y）y'＝（x－y） では、HOMOGENEOUS(（x－y）/（x+y）, x, y, x0, y0)として、
y = IF(x > 0, － √2√(x^2 + 1) － x) ∨ y = IF(x > 0, √2√(x^2 + 1) － x)が得られる。

（4） p(x, y) + q(x, y)y' = 0、ただし、dp/dy＝ dq/dx－－－ EXACT(p, q, x, y, x0, y0) または EXACT_GEN(p, q, x, y, c)

たとえば、(x+y^2)＋2xyy'＝0 のときは、dp/dy＝2y、dq/dx＝2yなので、EXACT(x + y^2, 2xy, x, y, x0, y0)と入力して、yについて解けば、
y = － √2√((x0(x0 + 2y0^2) － x^2)/x)/2 ∨ y = √2√((x0(x0 + 2y0^2) － x^2)/x)/2

（5）ベルヌイの微分方程式 BERNOULLI_ODE(p, q, k, x, y, x0, y0) またはBERNOULLI_ODE_GEN(p, q, k, x, y, c)

Bernoulli（ベルヌイ（ベルヌーイとも書く））は、ヨハン・ベルヌーイ。懸垂線の研究などで知られているスイスの数学者（1667－1748）である。一族に優れた数学・物理学者を輩出している。
ベルヌイの微分方程式とは、y' + p(x)y = q(x)y^k、ここで、kは、2以上の整数（0または1の時は変数分離で解ける）とする。
たとえば、y'＋y＝xy^3では、BERNOULLI_ODE_GEN(1, x, 3, x, y, c)とすると、1/y^2 = ce^(2x) + x + 1/2が得られる。

（6）クレローの微分方程式 CLAIRAUT(p, q, x, y, v, c)

アレクシス・クロード・ド・クレロー（Alexis Claude de Clairault）は、フランスの数学者（1713年－1765年）。（クレーローとも書くようである）
p(x*v － y) = q(v)、ここで、v＝y’である。
たとえば、y＝x*v＋v^2の場合は、解は、[cx － y + c^2 = 0, x + 2v = 0]のようにベクトルで与えられる。
ここで、前項は、一般解であり、後項は、特異解である。この場合は、dy/dx＝v＝－x/2から、y＝－x^2/4が得られる。
特異解は、一般解の包絡線になっている。

（7） r(x, y)y' + p(x)b(y) = q(x)、ただし、(db/dy)/r がyを含まない　
　ALMOST_LIN(r, b, p, q, x, y, x0, y0)　または　ALMOST_LIN_GEN (r, b, p, q, x, y, c)

たとえば、（2y+1）y'＋x^2（y^2＋y）＝x^2のとき、ALMOST_LIN_GEN(2y + 1, y^2 + y, x^2, x^2, x, y, c)と入力すると
y^2 + y = ce^(－ x^3/3) + 1が得られる。

（8）積分因数を使って完全微分方程式に変形できるもの　
　INTEGRATING_FACTOR(p, q, x, y, x0, y0) または INTEGRATING_FACTOR_GEN(p, q, x, y, c)

完全微分方程式というのは、X（x,y）dx＋Y（x,y）dy＝0と書いたときに、ある関数φ（x,y）の全微分に左辺がなっているものをいう。
積分因数というのは、そのままでは、全微分になっていないがある関数を掛けた場合に全微分になるとき、その係数を積分因数という。
DERIVEの記法では、p(x, y) + q(x, y)y' = 0としているので、
たとえば、（x－y）dx＋xdy＝0は、（x－y）＋x*ｙ’＝0なので、INTEGRATING_FACTOR(x－y, x, x, y, x0, y0)とすると、
y = － xLN(x) + xLN(x0) + y0x/x0が得られる。

（9）単項テスト　MONOMIAL_TEST(p, q, x, y)

p(x, y) + q(x, y)y'＝0で、積分因数がx^m*y^nの形であるかどうかを調べる関数
（8）の例では、1/x^2が積分因数として得られる。このテストで因数が得られる場合は、それを掛けた式は、（4）のケースに帰着できる。

(10) y' = r(x, y)、ただし、t＝x/yとしたとき、td(t*r）/dx／d(t*r)/dyがx及びyを含まない
GEN_HOM(r, x, y, x0, y0) または GEN_HOM_GEN(r, x, y, c)

たとえば、ｙ’＝（（x/y）＋1）^2のようなもの。
この解は、複雑なので、ここでは、省略する。

(11) y'＝r（u）、ただし、u = px + qy + k、p、q、kは、x,yを含まない
　FUN_LIN_CCF(r, p, q, k, x, y, x0, y0)　または FUN_LIN_CCF_GEN(r, p, q, k, x, y, c)

たとえば、ｙ’＝（ax+by+c）、FUN_LIN_CCF(ax + by + c, a, b, c, x, y, x0, y0)とすると
y = e^(bx － bx0)(a(bx0 + 1) + b(by0 + c))/b^2 － (abx + a + bc)/b^2と解が求められる。

(12) y'＝ r((ax+by+d)/(px+qy+k))－－LIN_FRAC(r, a, b, d, p, q, k, x, y, x0, y0) または LIN_FRAC_GEN(r, a, b, d, p, q, k, x, y, c)

これは、(11)を含んでいる。係数a、b等は、定数である。
たとえば、（2x－y+3)－（x－2y+3)ｙ’＝0では、LIN_FRAC_GEN((2x － y + 3)/(x － 2y + 3), 2, －1, 3, 1, －2, 3, x, y, c)とすると
LN((x + 1)^2/(x^2 + x(3 － y) + y^2 － 3y + 3)) = 2LN(x + 1) + 2cとなるが、両辺のexpを取って整理すると、
x^2 － xy + 3x + y^2 － 3y = βとなる。βは、定数。（※（11）のFUN_LIN_CCF、FUN_LIN_CCF_GENを含んでいませんでした。2014/8/19）

(13) 擬線形方程式　r(x, y)y' + p(x)b(y) = q(x)、ただし、（db/dy）／rは、yを含まない
　ALMOST_LIN(r, b, p, q, x, y, x0, y0) 　または ALMOST_LIN_GEN (r, b, p, q, x, y, c)

たとえば、y^2ｙ’/x－y^3=0では、ALMOST_LIN(y^2/x, y^3, －1, 0, x, y, x0, y0)とすると、
実数解は、y = y0e^(x^2/2 － x0^2/2)が得られる。（※説明が（7）と重複していました。お詫びして訂正します。2014/8/19）