常微分方程数值方法1、欧拉方法：1,,1,0),,(1-=+=+n k y t hf y y k k k k .function E=euler(f,a,b,ya,n)% Input - f is the function entered as a string 'f'% - a and b are the left and right end points% - ya is the initial condition y(a)% - n is the number of steps% Output - E=[T' Y'] where T is the vector of abscissas and% Y is the vector of ordinatesh=(b-a)/n;T=zeros(1,n+1);Y=zeros(1,n+1);T=a:h:b;Y(1)=ya;for j=1:nY(j+1)=Y(j)+h*feval(f,T(j),Y(j));endE=[T' Y'];【例】 用欧拉方法求解区间]3,0[内的初值问题：1)0(,2'=-=y yt y 。

f=inline('(t-y)/2','t','y');a=0;b=3;ya=1;n=12;%n=3,6,12,24,48,96...E=euler(f,a,b,ya,n),plot(E(:,1),E(:,2),'r*'),hold on符号解：y=dsolve('Dy=(t-y)/2','y(0)=1')h=(3-0)/12;t=0:h:3;y=eval(y);[t' y']用图比较数值解：(f 为ode 函数文件)ode45('f',[0,3],1)2、休恩(Huen)方法(即改进Euler 方法)：1,,1,0)],,(,(),([211-=+++=++n k y t hf y t f y t f hy y k k k k k k k kfunction H=heun(f,a,b,ya,n)% Input - f is the function entered as a string 'f'% - a and b are the left and right end points% - ya is the initial condition y(a)% - n is the number of steps% Output - H=[T' Y'] where T is the vector of abscissas and% Y is the vector of ordinatesh=(b-a)/n;T=zeros(1,n+1);Y=zeros(1,n+1);T=a:h:b;Y(1)=ya;for j=1:nk1=feval(f,T(j),Y(j));k2=feval(f,T(j+1),Y(j)+h*k1);Y(j+1)=Y(j)+(h/2)*(k1+k2);endH=[T' Y'];【例】 用休恩方法求解区间]3,0[内的初值问题：1)0(,2'=-=y yt y 。

f=inline('(t-y)/2','t','y');a=0;b=3;ya=1;n=12;H=heun(f,a,b,ya,n),plot(H(:,1),H(:,2),'r*'),hold on用图比较数值解：(f 为ode 函数文件)ode45('f',[0,3],1)3、四阶龙格－库塔方法：)22(643211k k k k hy y k k ++++=+function R=rk4(f,a,b,ya,n)% Input - f is the function entered as a string 'f'% - a and b are the left and right end points% - ya is the initial condition y(a)% - n is the number of steps% Output - R=[T' Y'] where T is the vector of abscissas and% Y is the vector of ordinatesh=(b-a)/n;T=zeros(1,n+1);Y=zeros(1,n+1);T=a:h:b;Y(1)=ya;for j=1:nk1=h*feval(f,T(j),Y(j));k2=h*feval(f,T(j)+h/2,Y(j)+k1/2);k3=h*feval(f,T(j)+h/2,Y(j)+k2/2);k4=h*feval(f,T(j)+h,Y(j)+k3);Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6;endR=[T' Y'];【例】 用4阶龙格－库塔方法求解区间]3,0[内的初值问题：1)0(,2'=-=y yt y 。

f=inline('(t-y)/2','t','y');a=0;b=3;ya=1;n=3;R=rk4(f,a,b,ya,n)plot(R(:,1),R(:,2),'r*'),hold on用图比较数值解：(f 为ode 函数文件)ode45('f',[0,3],1)4、阿当姆斯－巴什弗斯－摩尔顿方法：（预测－校正公式）⎪⎪⎩⎪⎪⎨⎧++-+=+-+-+=+--+---)),(9195(24)5559379(241121123p x f f f f h y y ff f f h y p k k k k k k k k k k kfunction A=adams(f,T,Y)% Input - f is the function entered as a string 'f'% - T is the vector of abscissas% - Y is the vector of ordinates% Remark:The first four coordinates of T and% Y must have starting values obtained with RK4% Output - A=[T' Y'] where T is the vector of abscissas and% Y is the vector of ordinatesn=length(T);if n<5,break,end;F=zeros(1,4);F=feval(f,T(1:4),Y(1:4));h=T(2)-T(1);for k=4:n-1% Predictorp=Y(k)+(h/24)*(F*[-9 37 -59 55]');T(k+1)=T(1)+h*k;F=[F(2) F(3) F(4) feval(f,T(k+1),p)];% CorrectorY(k+1)=Y(k)+(h/24)*(F*[1 -5 19 9]');F(4)=feval(f,T(k+1),Y(k+1));endA=[T' Y'];【例】 求解区间]3,0[内的初值问题：1)0(,2'=-=y y t y 。

f=inline('(t-y)/2','t','y'); T=zeros(1,25);Y=zeros(1,25);T=0:1/8:3; Y(1:4)=[1 0.94323919 0.89749071 0.86208736];A=adams(f,T,Y)5、米尔尼－辛普生（Milne-simpson ）方法：预测:)22(341231k k k k k f f f h y p +-+=---+ 修正：292811k k k k p y p m -+=++，),(111+++=k k k m t f f 校正:)4(31111+--++++=k k k k k f f f h y y function M=milne(f,T,Y)% Input - f is the function entered as a string 'f'% - T is the vector of abscissas% - Y is the vector of ordinates% Remark:The first four coordinates of T and% Y must have starting values obtained with RK4% Output - M=[T' Y'] where T is the vector of abscissas and% Y is the vector of ordinatesn=length(T);if n<5,break,end;F=zeros(1,4);F=feval(f,T(1:4),Y(1:4));h=T(2)-T(1);pold=0;yold=0;for k=4:n-1% Predictorpnew=Y(k-3)+(4*h/3)*(F(2:4)*[2 -1 2]');% Modifierpmod=pnew+28*(yold-pold)/29;T(k+1)=T(1)+h*k;F=[F(2) F(3) F(4) feval(f,T(k+1),pmod)];% CorrectorY(k+1)=Y(k-1)+(h/3)*(F(2:4)*[1 4 1]');pold=pnew;yold=Y(k+1);F(4)=feval(f,T(k+1),Y(k+1));endM=[T' Y'];【例】 求解区间]3,0[内的初值问题：1)0(,2'=-=y y t y 。

f=inline('(t-y)/2','t','y'); T=zeros(1,25);Y=zeros(1,25);T=0:1/8:3; Y(1:4)=[1 0.94323919 0.89749071 0.86208736];M=milne(f,T,Y)6、哈明(Hamming)方法：预测：)22(341231k k k k k f f f h y p +-+=---+ 校正：)),(2(83)9(8111121++--+++-++-=k k k k k k k p t f f f h y y y function H=hamming(f,T,Y)% Input - f is the function entered as a string 'f'% - T is the vector of abscissas% - Y is the vector of ordinates% Remark:The first four coordinates of T and% Y must have starting values obtained with RK4% Output - H=[T' Y'] where T is the vector of abscissas and% Y is the vector of ordinatesn=length(T);if n<5,break,end;F=zeros(1,4);F=feval(f,T(1:4),Y(1:4));h=T(2)-T(1);pold=0;cold=0;for k=4:n-1% Predictorpnew=Y(k-3)+(4*h/3)*(F(2:4)*[2 -1 2]');% Modifierpmod=pnew+112*(cold-pold)/121;T(k+1)=T(1)+h*k;F=[F(2) F(3) F(4) feval(f,T(k+1),pmod)];% Correctorcnew=(9*Y(k)-Y(k-2)+3*h*(F(2:4)*[-1 2 1]'))/8;Y(k+1)=Y(k-1)+(h/3)*(F(2:4)*[1 4 1]');pold=pnew;cold=cnew;F(4)=feval(f,T(k+1),Y(k+1));endH=[T' Y'];【例】 求解区间]3,0[内的初值问题：1)0(,2'=-=y y t y 。