2、,y(i-1));k2=feval(ode,t(i),y(i-1)+h*k1);y(i)=y(i-1)+h*(k1+k2)/2;end三龙格—库塔法(Runge-Kutta)欧拉公式可改写为它每一步计算f(xi,yi)一次，截断误差为O(h2)标准四阶龙格—库塔公式每一步计算f(x,y)四次，截断误差为O(h5)01/21/21/201/210011/62/62/61/6对于两个分量的一阶常微分方程组用经典4阶Runge-Kutta法求解的格式为n级显式Runge-Kutta方法的一般计算格式：其中Adams外插公式（Adams-Bashforth公式）是一
3、类k+1步k+1阶显式方法三步法（k=2）,四步法（k=3）,Adams内插公式（Adams-Moulton公式）是一类k+1步k+2阶隐式方法三步法（k=2）,Adams预估-校正方法（Adams-Bashforth-Moulton公式）一般取四步外插法与三步内插法结合。#include#include#include#defineTRUE1main(){intnstep_pr,j,k;floath,hh,k1,k2,k3,k4,t_old,t_limit,t_mid,t_new,t_pr,y,ya,y
4、n;doublefun();printf("Fourth-OrderRunge-KuttaScheme");while(TRUE){printf("Intervaloftforprinting?");scanf("%f",&t_pr);printf("Numberofstepsinoneprintinginterval?");scanf("%d",&nstep_pr);printf("Maximumt?");scanf("%f",&t_limit);y=1.0;/*Settingtheinitialvalueofthesolution*/
5、h=t_pr/nstep_pr;printf("h=%g",h);t_new=0;/*Timeisinitialized.*/hh=h/2;printf("--------------------------------------");printf("ty");printf("--------------------------------------");printf("%12.5f%15.6e",t_new,y);do{for(j=1;j<=nstep_pr;j++){t_old=t_new;t_new=t_new+h;yn=y;t
6、_mid=t_old+hh;yn=y;k1=h*fun(yn,t_old);ya=yn+k1/2;k2=h*fun(ya,t_mid);ya=yn+k2/2;k3=h*fun(ya,t_mid);ya=yn+k3;k4=h*fun(ya,t_new);y=yn+(k1+k2*2+k3*2+k4)/6;}printf("%12.5f%15.6e",t_new,y);}while(t_new<=t_limit);printf("--------------------------------------");printf("Maximumtlimitex
7、ceeded");printf("Type1tocontinue,or0tostop.");scanf("%d",&k);if(k!=1)exit(0);}}doublefun(y,t)floaty,t;{floatfun_v;fun_v=-y;return(fun_v);}四误差的控制我们常用事后估计法来估计误差，即从xi出发，用两种办法计算y(xi+1)的近似值。记为从xi出发以h为步长得到的y(xi+1)的近似值，记为从xi出发以h/2为步长跨两步得到的y(xi+1)的近似值。设给定精度为ε。如果不等式成立，则即为y(xi+1)的满足精度要求的
8、近似值。自适应：使用2个不同的h。如果一个大的h和一
