maximaの初歩的な操作18･･･１階常微分方程式

== maximaの初歩的な操作18 ==
･･･１階常微分方程式･･･
○Xmaxima，wxMaxima のインストール方法，基本操作については［この頁］参照
○wxmaximaを起動すると，初めは入力行も何もないので，とりあえず「何かキーを押します：例えばEnterキーを押す」と入力行が現れます．
○wxmaximaでは，入力行の末尾でShift+Enterを押します．

○wxMaximaで１階または２階の常微分方程式の一般解を求めるには，

(1)　ode2( )関数を使う方法
…ordinary differential equation of 1st and 2nd order
(2)　ラプラス変換を使ってdesolve( )関数で解く方法

があります．この頁では(1)のode2( )関数を使う方法を扱います．

１階微分方程式では積分定数は%cで表され，
２階微分方程式では積分定数は%k1 , %k2で表されます．

○wxmaximaで，

１階微分方程式で初期条件に対応する特殊解を求めるにはic1( )関数が，
…initial value problems for 1st order differential equations
２階微分方程式で初期条件に対応する特殊解を求めるにはic2( )関数が
…initial value problems for 2nd order differential equations
使えます．

※微分方程式の厳密解は求め方が知られていない場合や，解けても特殊関数の複雑な組合せとなる場合があります．
wxmaximaのヘルプ情報（英文）によれば，
１階常微分方程式については，

線形，変数分離形，同次形，完全微分形，ベルヌーイ形などが試され，解けなければfalseが返される

とされています．
２階常微分方程式については，

定数係数，定数係数に変換できる同次線形，オイラーの微分方程式，定数変化法によって解けるもの，２つの１階微分方程式に分けられるものなどが試されます．解けなければfalseが返される

　級数を用いた解については，このような制限はあまりないがこの頁では扱いません．

≪頭の片隅に！≫
　厳密解･･･有限回の代数的変形や積分操作によって得られる解･･･は，コンピュータでも求められない場合がある．（変数変換のアイデア次第で解ける問題は，人の方が強い場合がある）

　解き方の分かっている問題に当てめる場合や，理論上無限回（実際上は巨大な回数）の操作を用いる級数解や数値解を求める場合コンピュータの方が強い．

■１階常微分方程式

≪例１≫
【微分方程式】
$\frac{dy}{dx}=x^2y^2$
【wxMaximaでの入力】
ode2( 'diff(y , x)=x^2*y^2 , y , x);

【出力される結果】
$-\frac{1}{y}=\frac{x^3}{3}+ %c$ …(*)

【補足説明】
○関数ode2( )の引数は
ode2( 微分方程式 , 従属変数 , 独立変数);
の順に指定します．

メニューから対話型メニューで行うには，「方程式」→「微分方程式を解く」と進んで，上記の引数３個を書き込みます．

○その場合において，注意すべきこととして導関数を表す記号diff(y , x)にはシングルクォートを付けて'diff(y , x)と書きます．（単にdiff(y , x)と書くと，この文脈の中では文字yを文字xで（偏）微分した0に変わります．ここでは，この式自体を生かしたいので，名詞化処理処理としてシングルクォートを付けるようです．）
　次のように，xの関数であることを明示して使う場合には，シングルクォートを付けても付けなくても同じ結果が得られます．
ode2( diff(f(x),x)=x^2*f(x)^2 , f(x) , x);

○微分方程式の部分で等号を省略して単に式を指定した場合は，=0が省略されたものとみなされます．（右辺＝0）

○diff( )の第３引数の1を省略すれば，第１次導関数を表しますが，第２次導関数を表す場合はdiff(y , x , 2 )のように第３引数を2にします．

○この一般解は，次のようにyについて解かれた形に直せます．
$y=-\frac{3}{x^3+ 3C}$

○一般解(*)のうちで，与えられた初期条件「x=0のときy=−3」を満たす特殊解を求めるには，
ic1( 一般解 , x=0 , y=−3);
とします．
$-\frac{1}{y}=\frac{x^3+ 1}{3}$ になります．

　このように，先に一般解を求めてから，続いてある初期条件を満たす特殊解を求めるには，一般解の式に名前を付けておいて後でそれを引用するとよいので，上記の一般解を求めるコマンドを
siki1 : ode2( 'diff(y , x)=x^2*y^2 , y , x);
などとしておいてから
ic1( siki1 , x=0 , y=−3);
とすればよいことになります．
　Maximaでは単に%と書けば直前の出力式を表すことになっているので，一般解を求めた直後に初期条件を加える場合には，式にな名前がなくても
ic1( % , x=0 , y=−3);
でできます．

メニューから対話型メニューで行うには，「方程式」→「特殊解を求める（1階微分方程式）」と進んで，上記の引数３個を書き込みます．

≪例２≫
【微分方程式】
$x\frac{dy}{dx}+ 2y=e^x$
$x=1$ のとき $y=0$
【wxMaximaでの入力】
siki2 : ode2(x*'diff(y,x)+2*y=%e^x , y , x);
ic1(siki2,x=1,y=0);

【出力される結果】
$y=\frac{(x-1)e^x+ C}{x^2}$

Maximaでは積分定数は%cで表わし，自然対数は%eで表します

$y=\frac{(x-1)e^x}{x^2}$

≪例３≫
【微分方程式】
$(1-x^2)\frac{dy}{dx}+ xy=2x$
$x=0$ のとき $y=1$
【wxMaximaでの入力】
不定積分の結果として対数が登場するときに，その真数に絶対値を付けるには，あらかじめ変数logabsのプロパティをtrueにしておきます
logabs : true; …(1)
siki3 : ode2((1-x^2)*'diff(y,x)+x*y=2*x, y, x); …(2)
累乗根（radicals），指数関数，対数関数が含まれる式を標準的な（canonical ）式に直すには，関数radcan( )を実行します．
siki4 : radcan(siki3); …(3)
ic1(siki4 , x=0 , y=1); …(4)

【補足説明】
○この微分方程式の解が対数関数を含む式になるということは，やってから分かることで，やってみないと分かりません．だから，(2)よりも先に(1)を実行しておくというのは実際には無理です．まず(2)を実行して対数関数が登場するのを見届けたら，(1)を実行してからもう一度(2)をやり直すことになります．

　この煩わしさを避けるには，
siki3 : ode2((1-x^2)*'diff(y,x)+x*y=2*x, y, x) , logabs; …(2’)
　のようにカンマで区切ってlogabsを続けて実行するとよいでしょう．（この場合，logabsの値はこのコマンドを実行するとき trueとなり，以後の設定に影響しません･･･とはいえ，logが登場するかどうかが，やってみないと分からなければ，少し違うだけとも言えますが）

○(3)について
　例えば $\log e^x=x$ ， $e^{\log x}=x$ のように累乗根，指数関数，対数関数の組合せによってはもっと簡単な表現式がある場合に，式を簡単にする関数です．
　この関数は，メニューから「式の変形」→「関数の整理」と進んでも実行することができます．

【出力される結果】
(2)→ $y=\Bigl(2e^{-\frac{log(|x^2-1|)}{2}} + C\Bigr)e^{\frac{log(|x^2-1|)}{2}}$

Maximaでは積分定数は%cで表わし，自然対数は%eで表します

(3)→ $y=C\sqrt{|x^2-1|}+ 2$
(4)→ $y=2-\sqrt{|x^2-1|}$

≪例４≫
【微分方程式】
$(x+ y)^2\frac{dy}{dx}=1$
【wxMaximaでの入力】
ode2((x+y)^2*'diff(y,x)=1 , y , x); …(*)

【出力される結果】
false

【補足説明】
　未知関数yやその導関数y’の積などが登場する非線形微分方程式をそのままコンピュータにかけても厳密解が求まらないことがあります．このような場合でも，何らかの変数変換を行えば解けることがあります．

⇒　［コンピュータが働くように変形する］
$x+ y=t$ とおくと
$y=t$ $-x\rightarrow$ $\frac{dy}{dx}=\frac{dt}{dx}-1$
$t^2(\frac{dt}{dx}-1)=1$
$\frac{dt}{dx}-1=\frac{1}{t^2}$
$\frac{dt}{dx}=\frac{1}{t^2}+ 1$ は変数分離形→解ける
【wxMaximaでの入力】
ode2(t^2*('diff(t,x)-1)=1, t, x); …(**)
【出力される結果】
$t-\tan^{-1}t=x+ C$
変数を元に戻すと
$x+ y-\tan^{-1}(x + y)=x+ C$
$y=\tan^{-1}(x + y)+ C$ …（答）

その他：１階常微分方程式の例と結果

○≪変数分離形は解ける…非線形でも可≫

○微分方程式

wxMaximaでの入力

出力される結果（および備考）

○ $\frac{dy}{dx}=\frac{xy}{1+ 2x^2}$

ode2( 'diff(y,x)=x*y/(1+2*x^2) , y , x);

$y=Ce^{\frac{\log(2x^2 + 1)}{4}}$
radcan(%);を実行することにより
$y=C(2x^2 + 1)^{\frac{1}{4}}$

○ $\frac{dy}{dx}-e^x y=2e^{2x}$

ode2( 'diff(y,x)-%e^x*y=2*%e^(2*x) , y , x);

$y=(2(-e^x-1)e^{-e^{x}}+ C)e^{e^x}$
radcan(%);を実行することにより
$y=Ce^e^x-2e^x-2$

○ $\frac{dy}{dx}=\frac{(x-a)(y-p)}{(x-b)(y-q)}$

ode2( 'diff(y,x)=(x-a)*(y-p)/((x-b)*(y-q)) , y , x);

$(p-q)\log(y-p)+ y=(b-a)\log(x-b)+ x+ C$

○≪同次形（ $\frac{dy}{dx}=f(\frac{y}{x})$ の形になるもの）は解ける…非線形でも可≫

○微分方程式

wxMaximaでの入力

出力される結果（および備考）

○ $x\frac{dy}{dx}=x+ y$

ode2(x*'diff(y,x)=x+y , y , x);

$y=x(\log x + C)$
この解には対数関数が登場するので，既に述べた内容を踏まえると，あらかじめlogabs : true;を実行しておくべしとなるが，なぜか？そうすると解け切らずに
$y=|x|\Bigl(\int \frac{1}{|x|}dx+ C \bigr)$
などという結果が表示されます．
　期待される答は， $y=x(\log|x|+ C)$ であるはずだが?

○ $x^2+ y^2=2xy\frac{dy}{dx}$

ode2(x^2+y^2=2*x*y*'diff(y,x) , y , x);

$-\frac{y^2-x^2}{x}=C$

○ $(2xy-x^2)\frac{dy}{dx}+y^2-2xy=0$

ode2((2*x*y-x^2)*'diff(y,x)+y^2-2*x*y=0 , y , x);

$xy^2-x^2y=C$

○≪線形微分方程式（ $y'+ P(x)=Q(x)$ の形）は解ける≫

○微分方程式

wxMaximaでの入力

出力される結果（および備考）

○ $\frac{dy}{dx}-2y=e^{5x}$

ode2( 'diff(y,x)-2*y=%e^(5*x) , y , x);

$y=e^{2x}\Bigl( \frac{e^{3x}}{3}+ C\Bigr)$
式の展開expand(%);を実行すれば
$y=\frac{e^{5x}}{3} + Ce^{2x}$

○ $\frac{dy}{dx}+ y=cos x$

ode2( 'diff(y,x)+y=cos(x) , y , x);

$y=e^{-x} \Bigl( \frac{e^x(\sin x+ \cos x)+ C}{2}\Bigr)$
式の展開expand(%);を実行すれば
$y=\frac{\sin x}{2}+ \frac{\cos x}{2} + Ce^{-x}$

○ $(e^y-x)\frac{dy}{dx}=y$

ode2((%e^y-x)*'diff(y,x)=y , y , x);

$e^y-xy=C$
元の微分方程式は線形微分方程式とは言えないが，独立変数と従属変数を逆にすると
$\frac{dx}{dy}+ (\frac{1}{y})x=\frac{e^y}{y}$
の形になるので（少し計算すれば分かる），xがyの関数として線形微分方程式で表されることになる．このように立場を変えれば解ける問題は自動的に試されるようです．

○≪完全微分形（ $\frac{\partial F}{\partial x}dx+\frac{\partial F}{\partial y}dy=0$ の形）は解ける≫

○微分方程式

wxMaximaでの入力

出力される結果（および備考）

○ $(2x+ y)dx+xdy=0$

ode2(2*x+y+x*'diff(y,x)=0 , y, x);

$y=\frac{C-x^2}{x}$
微分形式dx ,dyのままode2()関数に記入するのは無理なようなので，
$(2x+ y)+ \frac{dy}{dx}=0$
の形で方程式を解きます．
　この方程式は $F(x,y)=x^2+ xy$ とおくと
$\frac{\partial F}{\partial x}=2x+ y,\frac{\partial F}{\partial y}=x$
が成り立つので完全微分形になっており， $F(x,y)=x^2+ xy=C$ が一般解です．wxMaximaから得られた解を少し変形すればこの式になります．
　ただし，この問題は同次形ともいえるのでMaximaが同次形を解く手順で解いている可能性があります．

○ $ye^xdx+(e^x-\cos y)dy=0$

ode2(y*%e^x+(%e^x-cos(y))*'diff(y,x)=0 , y , x);

$e^xy-\sin y=C$
微分形式dx ,dyのままode2()関数に記入するのは無理なようなので，
$ye^x+(e^x-\cos y)\frac{dy}{dx}=0$
の形で方程式を解きます．
この方程式は $F(x,y)=e^xy-\sin y$ とおくと
$\frac{\partial F}{\partial x}=e^xy,\frac{\partial F}{\partial y}=e^x- \cos y$
が成り立つので完全微分形になっており， $F(x,y)=e^xy-\sin x=C$ が一般解です．

○≪全微分方程式で積分因子の発見を要する問題は，解ける場合と解けない場合がある≫

○微分方程式

wxMaximaでの入力

出力される結果（および備考）

○ $(x^2 + 2x-y^2)dx- 2ydy=0$

ode2(x^2+2*x-y^2-2*y*'diff(y,x)=0 , y , x);

$x^2e^x-e^xy^2=C$
　この方程式は $F(x,y)=e^x(x^2-y^2)$ とおくと
$\frac{\partial F}{\partial x}=e^x(x^2-y^2)+e^x(2x)=e^x(x^2+ 2x-y^2)$
$\frac{\partial F}{\partial y}=e^x(-2y)$
が成り立つので積分因子 $M(x,y)=e^x$ を掛ければ完全微分形になります．
この問題では，積分因子M(x,y)がxだけの関数である場合に，Maximaがこれを見つけて解けます．

○ $2xydx+(2x^2-3y)dy=0$

ode2(y*%e^x+(%e^x-cos(y))*'diff(y,x)=0 , y , x);

$e^xy-\sin y=C$
　この方程式は $F(x,y)=x^2y^2-y^3$ とおくと
$\frac{\partial F}{\partial x}=2xy^2=y(2xy)$
$\frac{\partial F}{\partial y}=2x^2y-3y^2=y(2x^2-3y)$
が成り立つので積分因子 $M(x,y)=y$ を掛ければ完全微分形になります．
この問題では，積分因子M(x,y)がyだけの関数である場合に，Maximaがこれを見つけて解けます．

○ $(x^2y^2-y)dx+(x^2y^2-y)dy=0$

ode2(x^2*y^2-y+(x^2*y^2-x)*'diff(y,x)=0 , y , x);

$false$
　この方程式は $F(x,y)=x+y+\frac{1}{xy}$ とおくと
$\frac{\partial F}{\partial x}=1-\frac{1}{x^2y}=\frac{1}{x^2y^2}(x^2y^2-y)$
$\frac{\partial F}{\partial y}=1-\frac{1}{xy^2}=\frac{1}{x^2y^2}(x^2y^2-x)$
が成り立つので積分因子 $M(x,y)=\frac{1}{x^2y^2}$ を掛ければ完全微分形になります．
しかし，この問題では積分因子M(x,y)がxyの関数である場合に，Maximaがこれを見つけて解くことはできません．

ode2(1/(x^2*y^2)*(x^2*y^2-y)+1/(x^2*y^2)*(x^2*y^2-x)*'diff(y,x)=0
, y , x);

の形に変形すれば，次の解が得られます．

$\frac{xy^2+ x^2y+ 1}{xy}=C$

○ $y(xy+ y)dx+ x(6xy-4)dy=0$

ode2(y*(x*y+1)+x*(6*x*y-4)*'diff(y,x)=0 , y , x);

$x=Ce^{\frac{2\log(xy-1)+ 4\log(xy)}{5}}$

radcan(%);を実行することにより

$x=Cx^{\frac{4}{5}}y^{\frac{4}{5}}(xy-1)^{\frac{2}{5}}$

　この解は少し変形すると

$x^{\frac{1}{2}}y^3-x^{-\frac{1}{2}}y^2=C$

と書けます．

　この方程式を全微分方程式として解くには，積分因子 $M(x,y)=\frac{y}{x\sqrt{x}}$ を見つける必要がありますが，この問題では積分因子M(x,y)がx,yの関数である場合でも，Maximaがこれを見つけて解くことができます．

○≪ベルヌーイ形（ $y'+ P(x)y=Q(x)y^n(n\neq 0,1$ の形）は解ける≫

○微分方程式

wxMaximaでの入力

出力される結果（および備考）

○ $\frac{dy}{dx}+ \frac{y}{2}=y^3\cos x$

ode2( 'diff(y,x)+y/2=y^3*cos(x) , y , x);

$y=\frac{e^{-\frac{x}{2}}}{C-e^{-x}(\sin x-\cos x)}$
この解は変形して次の形に書くことができます
$y^2=\frac{1}{\cos x-\sin x + Ce^x}$

○ $2xy^2\frac{dy}{dx}-y^3+ x=0$

ode2(2*x*y^2*'diff(y,x)-y^3+x=0 , y , x);

$\frac{2y^3-6x}{3x^{\frac{3}{2}}}=C$

○≪変数変換によって既知の形に直す問題は，解ける場合と解けない場合がある≫

○微分方程式

wxMaximaでの入力

出力される結果（および備考）

○ $\frac{dy}{dx}=\frac{4x-2y+ 3}{2x-y+ 1}$

ode2((2*x-y+1)*'diff(y,x)=4*x-2*y+3 , y , x);

$\frac{y^2 +(-4x-2)y+ 4x^2+ 6x}{2}=C$
この問題は $2x-y+ 1=z$ とおくと，変数分離形になり
$\frac{dz}{dx}=-\frac{1}{z}\rightarrow \frac{z^2}{2}=-x+ C_1 \rightarrow z^2=-2x+ C_2$
$\rightarrow z^2+ 2x=C_2 \rightarrow (2x-y+ 1)^2+ 2x=C$
ですが，maximaはこの変形を見つけるようです．

○ $\frac{dy}{dx}=\frac{1}{(x-y)^2}$

ode2((x-y)^2*'diff(y,x)=1 , y , x);

$false$
この問題は $x-y=z$ とおくと，変数分離形になり
$1-\frac{dz}{dx}=\frac{1}{z^2}\rightarrow \cdots \rightarrow \frac{x-y-1}{x-y+ 1}=Ce^{2y}$
ですが，maximaはこの変形を見つけられないようです．