微分方程式

微分方程式

東大の後期を除けば，入試にはあまり出ないのではないかと思われるが，範囲としては面白いし，グラフィックも面白いので実験しよう。
（常）微分方程式は Ordinary Didderential Equation（ODE）
パッケージは plotdf（plot direction fields）

(%i1) load(plotdf)$

微分方程式を解くには，フーリエ変換を利用した desolve という関数と，陰関数もできる ode2（1,2階の常微分方程式を解く）という２つの関数がある。
desolve の初期値は atvalue であらかじめ与えておく必要があり，ode2のほうは初期値問題として ic1（initialconstantかな）またはic2（2階の方程式の初期値）という関数で解く。

「接線の傾きがその点のx座標の2倍で(0,1)を通る関数」を求める。

(%i2) atvalue(f(x),x=0,1)$desolve(diff(f(x),x)=2*x,f(x));
Result

(%i4) ode2('diff(y,x)=2*x,y,x);
Result

’は評価するという関数，これがないと diff(y,x)=0

(%i5) ic1(%,x=0,y=1);
Result

微分方程式と初期値を与えて，解を求める関数を作った。sode1（１階の常微分方程式の特殊解）

(%i6) sode1(e,a,b):=ic1(ode2(e,y,x),x=a,y=b)$

(%i7) sode1('diff(y,x)=2*x,0,1);
Result

これをベクトル場として目で見ることができて，初期値はマウスで指定する。すると赤い軌道がえがかれる。

(%i8) plotdf(2*x,[trajectory_at,0,1])$

Figure 1:
Result

変化量がそのy座標の-1/2倍

(%i9) desolve(diff(f(x),x)=-f(x)/2,f(x));
Result

(%i10) ode2('diff(y,x)=-y/2,y,x);
Result

(%i11) sode1('diff(y,x)=-y/2,0,1);
Result

(%i12) plotdf(-y/2)$

Figure 2:
Result

陰関数っぽいやつは，desolve ではうまくいかない。

(%i13) desolve(diff(f(x),x)=-x/f(x),f(x));
Result

(%i14) ode2('diff(y,x)=-x/y,y,x);
Result

(%i15) sode1('diff(y,x)=-x/y,0,1);
Result

(%i16) plotdf(-x/y)$

Figure 3:
Result

入試問題集から選んでみたが，軽く解ける。

(%i17) sode1('diff(y,x)=(2*y-3)*x,0,1);
Result

(%i18) ode2((y+'diff(y,x))*sin(x)=y*cos(x),y,x);
Result

(%i19) sode1('diff(y,x)+2*y=2*x*%e^(-x),0,0);
Result

(%i20) sode1(sqrt(%e^x+1)*('diff(y,x)+y)=1,0,0);
Result

(%i21) ode2('diff(y,x,2)-3*'diff(y,x)-4*y,y,x);
Result

(%i22) ic2(%,x=0,y=1,'diff(y,x)=9);
Result

2階の常微分方程式と初期値（yとy'）を与えて特殊解を出す関数を作った。

(%i23) sode2(e,a,b,c):=(ic2(ode2(e,y,x),x=a,y=b,'diff(y,x)=c))$

(%i24) sode2('diff(y,x,2)-3*'diff(y,x)-4*y,0,1,9);
Result