SCILABで微分方程式を解く。


単純な微分方程式 y'=2x の解は y=x^2

となるわけで、これをSCILABで解いて見る。この「解く」というのは積分を求めて定性的に解くわけではなく、微分≒差分を利用して 解析的に解くわけである。
SCILABでこれを解くのはode (ordinary differential equation solver)関数である。

まず微分方程式を
y' = f(x,y)

の形に書く。ここで

y'=f(x)

とは書けず、必ず

y'=f(x,y)

の形にしないと何故かODEを通せない。
y'=f(x,y)=2x なので

ydot = df(x,y)

と宣言して、中身を
ydot = 2*x

と定義する。これをSCILABでは

deff("ydot=df(x,y)" , "ydot=2*x")

と書く。
このydotとかdfとかx,yの変数名、関数名は任意である。
ちなみにSCILABの「行の最後が..であれば、次の行に続く」という機能を利用して

deff( ..
"ydot=df(x,y)" , ..
"ydot=2*x" ..
)

と書くのもオツだ。(笑)

それから xの初期値、x0の時のyの値、y0=y(x)を準備しておく
x0=0 ; y0=0 ;
次に、求めたいxの値の配列を作り
x=[0:1:10]; //0〜10まで1ずつ増やす

deを使って
y=ode(x0,y0,x,df);
とすると、yにy(x=0) , y(x=1) , y(x=2) ... x(x=10)
の値が入る。

deff("ydot=df(x,y)" , "ydot=2*x")

x0=0 ; y0=0 ; x=0:1:10;

y=ode(y0 , x0 , x , df );

plot(t,y)


ここで注目すべきは、一般に微分方程式をこの様に数値解析で求める時には RungeKutta法を基本として、早い話しが微分を差分で近似して求めるので
刻み幅(ここでは1)を小さくしないと、解が正しくないのだが、 odeは(多分)内部で刻み幅を自動的に決定し、その刻み幅で解を求め たのち、与えられた刻み幅の時の値を返すので、
x=[0:0.0001:10]; //0〜10まで0.0001ずつ増やす
でも
x=[0:5:10]; //0〜10まで5ずつ増やす
でも与えた個所のX、例えばy(x=10)の値は両方とも殆ど変わらないということだ。

(この微分を差分で〜という個所については、MaTXの項で書いているので参考にしてください)

MaTXの項で書いた、振動の微分方程式をSCILABで解くと

//---------------------------
deff( .. //deff宣言 "xdot=df(t,x)" , .. x'=df't,x)
"M=[[0 , 1] ; .. Mの行列
[ -(2*%pi)^2 , 0]]; ..
xdot=M*x ; " .. xは [x1,x2]'の形
)
//------------------------------------------------
//1行で書くなら deff("xdot=df(t,x)" , "M=[[0, 1];[ -(2*%pi)^2 , 0]]; xdot=M*x ; ")

x0=[1,0]' ; t0=0 ; t=0:0.01:5;
y=ode(x0,t0,t,df);
// y(1,:)はyの1行目だけを取り出したベクトル plot2d(t,y(1,:),5,"111","ode",[0,-1,5,1],[2,5,5,2]);xgrid(2) // [0,-1,5,1]はグラフの Xmin , Ymin , Xmax , Ymax // [2,5,5,2 ]はグラフの X補助目盛数、X主目盛数、Y補助目盛数、Y主目盛数、


return