igorで調和振動子のルンゲ・クッタ法による数値計算

ルンゲクッタ法を使って調和振動子の数値計算をしてみるという試み。

内容としてはかなり良くある話なんだけど、使うのがigor のマクロというのが普通と違うところ？

内容としてはC言語で書くのとほぼ同じなのだが、igorで書けば計算結果をすぐに描画できるかなというところ。
ただし、めちゃめちゃ計算は遅いのでメリットはあるのかよくわからない。

ルンゲ・クッタ法の備忘録として・・・。

一次元調和振動子の問題は、振り子の振動、バネの振動なんかでよく大学の物理でも出てくる問題である。
変移に比例した復元力が働く時の運動である。
微分方程式としては以下の通り。

macro ion_motion1D()

variable dt = 1e-8
variable xini = 1e-6
variable omegax = 2*pi*1e+6
variable nsize = 1000
variable i=1

make/D/N=(nsize)/O tim
make/D/N=(nsize)/O xion
make/D/N=(nsize)/O vion

//Rumge-Kutta parameter
variable dx1
variable dx2
variable dx3
variable dx4
variable dvx1
variable dvx2
variable dvx3
variable dvx4
Silent
DelayUpdate

xion[0]= xini
vion[0]=0
tim = dt*p

do

　dvx1=dt*(-omegax*omegax)*(xion[i-1]);
　dx1=dt*vion[i-1];
    dvx2=dt*(-omegax*omegax)*(xion[i-1]+1.0/2.0*dx1);
    dx2=dt*(vion[i-1]+1.0/2.0*dvx1);
    dvx3=dt*(-omegax*omegax)*(xion[i-1]+1.0/2.0*dx2);
    dx3=dt*(vion[i-1]+1.0/2.0*dvx2);
    dvx4=dt*(-omegax*omegax)*(xion[i-1]+dx3);
    dx4=dt*(vion[i-1]+dvx3);
    xion[i]=xion[i-1]+1.0/6.0*(dx1+2.0*dx2+2.0*dx3+dx4);
    vion[i]=vion[i-1]+1.0/6.0*(dvx1+2.0*dvx2+2.0*dvx3+dvx4);
    
    i=i+1
    //print i
    
while(i<=nsize)