計算する - ryamadaのコンピュータ・数学メモ

こちらでMATLABを使っている
フリーで何かあるかと調べるとScilabと言うのがあるらしい
Scilabのダウンロードとか解説については
手順
- Scilabから、適したＯＳのそれをダウンロードしてクリックオンリーでインストール
- こちらの微分方程式云々のところを試してみる

- 以下のソースは、はてなのスーパーpre記法の"Scilab"用(">|scilab|"で書き始めて"||<"で書き終える)を使用
- MATLABの場合は">|matlab|"で書き始めて"<||"で書き終える
  - こちらを
- Scilabのツールバーには"アプリケーション"→"MATLABからScilabへ変換"というコマンドがあるので、MATLABで書いたものをScilab化できそうだったのでやってみたけれど、できたのかできていないのか、未確認

function xdot=myfunc4(t,x)
 xdot=zeros(2,1);
 a=0.5;
 xdot(1)=-a*x(1)-sin(x(2));
 xdot(2)=x(1);
endfunction

x=-9:0.5:9;
v=-5:0.5:5;

dx=zeros(length(x),length(v),2);

for nx=1:length(x)
 for nv=1:length(v)
 dx(nx ,nv ,:)=myfunc4(0,[v(nv) x(nx)]);
 end
end

clf(); champ(x, v, dx(:,:,2), dx(:,:,1))

常微分方程式を解こう
- $\frac{dy}{dt}=y'= f(t,y)$ という微分方程式
- $\frac{dy}{dt}=y'=y^2-y\times \sin(t)+\cos(t)$
- 初期値、 $y0=0,t0=0$
- 計算すべきtの値 $t$ は0から0.1刻みで $2\pi$ まで
- ode()関数を使って微分方程式を解いてtの値に対応するyの値を求めて
- プロットする

function ydot=f(t, y),ydot=y^2-y*sin(t)+cos(t),endfunction
y0=0;t0=0;t=0:0.1:2*%pi;
y=ode(y0,t0,t,f)
plot(t,y)

- 同じことをRでやると

# package をインストールする
install.packages("deSolve")
# packageを読み込む
library(deSolve)
# 微分方程式を表す関数を作る
# ode()関数に渡す微分方程式を表す関数は、変数tと従属変数yと係数のベクトルとを順番に渡す
# ode()関数はparsという関数の係数などを引数で取らないといけないらしいので、以下のydotと定義する関数には使っていない引数parsがある
ydot<-function(t,y,pars){
# 返り値はリストで返す
# dy/dt = y^2-y*sin(t)+cos(t) という微分方程式
	list(y^2-y*sin(t)+cos(t))
}

# yの初期値を決める
y<-0
# 計算する時刻ポイントを決める
times<-seq(from=0,to=2*pi,by=0.1)
# ode()関数に処理させる
out1<-ode(y = y, times = times, func = ydot)
plot(out1)
par(new=TRUE)