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)
の値が入る。

ここで注目すべきは、一般に微分方程式をこの様に数値解析で求める時には RungeKutta法を基本として、早い話しが微分を差分で近似して求めるので
刻み幅（ここでは１）を小さくしないと、解が正しくないのだが、 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]'の形
)
//------------------------------------------------
//１行で書くなら 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,:)はｙの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