一阶常微分方程的初值问题求解

1、一阶常微分方程的初值问题 1.欧拉方法欧拉公式：算法（欧拉方法）（1）给定整数N；（2）；（3）（4）对k=2,3,N+1做：（5）返回基于欧拉方法的Matlab程序如下：function x,y=odeEuler(f,y0,a,b,h)n=(b-a)/h;x=a:h:b;y(1)=y0;for i=1:n y(i+1)=y(i)+h*feval(f,x(i),y(i);end2.改进的欧拉公式：算法（改进欧拉方法）：(1) 给定整数N；(2) ；(3)(4) 对k=2,3,N+1做：(5) 返回基于改进欧拉方法的Matlab程序如下：function x,y=odeIEuler(f,a,b,

2、y0,h)n=(b-a)/h;x=a:h:b;y(1)=y0;for i=1:n k1=h*feval(f,x(i),y(i); k2=h*feval(f,x(i+1),y(i)+k1); y(i+1)=y(i)+0.5*(k1+k2);end3.一类应用广泛的高精度的显式单步法-龙格-库塔（Runge-Kutta）方法，简称R-K方法标准四阶四段龙格-库塔公式算法（标准四阶四段龙格-库塔方法）（1）给定整数N；（2）；（3）（4）对k=2,3,N+1做：（5）返回基于标准四阶四段龙格-库塔方法Matlab程序如下：function x,y=oderk4(f,a,b,h,y0)n=(b-a)/

3、h;x=a:h:b;y(1)=y0;for i=1:n; k1=h*feval(f,x(i),y(i); k2=h*feval(f,x(i)+h/2,y(i)+k1/2); k3=h*feval(f,x(i)+h/2,y(i)+k2/2); k4=h*feval(f,x(i)+h,y(i)+k3); y(i+1)=y(i)+1/6*(k1+2*k2+2*k3+k4);end其中：f为常微分方程的右端项，a,求解区间的左右端点，h为自变量的步长，y0微分方程的初值，x,y分别为计算完成时的自变量取值和对应点上的函数值。4.标准四阶四段龙格-库塔方法应用实例 function z= f(x,y)z

4、=y-2x/y;function x,y=oderk4(f,a,b,h,y0)n=(b-a)/h;x=a:h:b;y(1)=y0;for i=1:n; k1=h*feval(f,x(i),y(i); k2=h*feval(f,x(i)+h/2,y(i)+k1/2); k3=h*feval(f,x(i)+h/2,y(i)+k2/2); k4=h*feval(f,x(i)+h,y(i)+k3); y(i+1)=y(i)+1/6*(k1+2*k2+2*k3+k4);end运行结果：>>x,y=oderk4(f,0,1,0.1,1)x = Columns 1 through 9 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 Columns 10 through 11 0.9000 1.0000y = Columns 1 through 9 1.0000 1.0954 1.1832 1.2649