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

一阶常微分方程的初值问题 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)); end 2.改进的欧拉公式： 算法（改进欧拉方法）： 给定整数N； ； 对k=2,3,…N+1做： 返回 基于改进欧拉方法的Matlab程序如下： function [x,y]=odeIEuler(f,a,b,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); end 3.一类应用广泛的高精度的显式单步法-----龙格-库塔（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)/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,求解区间的左右端点，hy0微分方程的初值，x,y分别为计算完成时的自变量取值和对应点上的函数值。 4.标准四阶四段龙格-库塔方法应用实例 function z= f(x,y) z=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.0000 y = Columns 1 through 9 1.0000 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 1.6125 Columns 10 through 11 1.6733 1.7321