09/07/14 13:01:33
>>56の続きです。
/* ルンゲクッタ法を1ステップ行う関数 */
/* x(t)=xin から x(t+h)=xout を求める */
void rungekutta(double xin[], double xout[], double t, double h)
{
int i;
double th2,h2,h6;
double dxdt0[N], dxdt1[N], dxdt2[N], dxdt3[N], xt[N];
double deriv(double, double [], double []);
h2=h*0.5;
h6=h/6.0;
th2=t+h2;
derivs(t,xin,dxdt0);
for(i=0;i<N;i++) xt[i] = xin[i] + h2*dxdt0[i];
derivs(th2,xt,dxdt1);
for(i=0;i<N;i++) xt[i] = xin[i] + h2*dxdt1[i];
derivs(th2,xt,dxdt2);
for(i=0;i<N;i++) xt[i] = xin[i] + h*dxdt2[i];
derivs(t+h,xt,dxdt3);
for(i=0;i<N;i++)
xout[i] = xin[i] + h6*(dxdt0[i]+2.0*(dxdt1[i]+dxdt2[i])+dxdt3[i]);
}