但是感觉依然不理想
|
/************************************************************************ * 用四阶(定步长)龙格--库塔法求解初值问题，其中一阶微分方程未y'=f(x,y) * 初始条件为x=x[0]时,y＝y[0]. * 输入: f--函数f(x,y)的指针 * x--自变量离散值数组(其中x[0]为初始条件) * y--对应于自变量离散值的函数值数组(其中y[0]为初始条件) * h--计算步长 * n--步数 * 输出: x为说求解的自变量离散值数组 * y为所求解对应于自变量离散值的函数值数组 ************************************************************************/ double runge_kuta(double(*f)(double,double),double x[], double y[],double h,int n) { int i; double xs,ys,xp,yp,dy; xs=x[0]+n*h; for(i=0;i<n;i++) { ys=y; dy=(*f)(x,y); //k1 y[i+1]=y+h*dy/6; xp=x+h/2; yp=ys+h*dy/2; dy=(*f)(xp,yp); //k2 y[i+1]+=h*dy/3; yp=ys+h*dy/2; dy=(*f)(xp,yp); //k3 y[i+1]+=h*dy/3; xp+=h/2; yp=ys+h*dy; dy=(*f)(xp,yp); //k4 y[i+1]+=h*dy/6; x[i+1]=xp; if(x[i+1]>=xs) return (0); } return(0); } |
