常微分方程的差分方法实验.docVIP

常微分方程的差分方法实验  实验目的及要求： （1）方法的，学会用解决某些实际问题，比较方法解题的不同处。 （2）熟悉Matlab编程环境，利用Matlab实现具体的。 用Matlab软件实现欧拉方法改进的欧拉方法龙格库塔方法，并用实例在计算机上计算。实验题目 （1）取h=0.1，用欧拉方法改进的欧拉方法四阶龙格-库塔方法求解初值问题： 并与精确解比较结果. ，取，计算。 实验过程： (1)f2 .m: function z=f2(x,y) z=1/(1+x^2)-2*(y^2); solvef2.m: function y=solvef2(x) y=x./(1+x.^2); Euler方法： 首先编写文件，如下： function E=Euler(f,a,b,N,ya) h=(b-a)/N; y=zeros(1,N+1); x=zeros(1,N+1); y(1)=ya; x=a:h:b; for i=1:N x(i+1)=x(i)+h; y(i+1)=y(i)+h*feval(f,x(i),y(i)); end E=[x,y]; 改进的Euler方法： 编写文件文件，如下： function E=MendEuler(f,a,b,N,ya) h=(b-a)/N; y=zeros(1,N+1); x=zeros(1,N+1); y(1)=ya; x=a:h:b; for i=1:N y1=y(i)+h*feval(f,x(i),y(i)); y2=y(i)+h*feval(f,x(i+1),y1); y(i+1)=(y1+y2)/2; end E=[x,y]; 四阶Runge-Kutta方法： 编写文件， function R=Rungkuta4(f,a,b,N,ya) h=(b-a)/N; x=zeros(1,N+1); y=zeros(1,N+1); x=a:h:b; y(1)=ya; for i=1:N k1=feval(f,x(i),y(i)); k2=feval(f,x(i)+h/2,y(i)+(h/2)*k1); k3=feval(f,x(i)+h/2,y(i)+(h/2)*k2); k4=feval(f,x(i)+h,y(i)+h*k3); y(i+1)=y(i)+(h/6)*(k1+2*k2+2*k3+k4); end R=[x,y]; (2)二阶Adams预报校正系统： function A=Adams2PC(f,a,b,N,ya) h=(b-a)/N; x=zeros(1,N+1); y=zeros(1,N+1); x=a:h:b; y(1)=ya; for i=1:N if i==1 y1=y(i)+h*feval(f,x(i),y(i)); y2=y(i)+h*feval(f,x(i+1),y1); y(i+1)=(y1+y2)/2; dy1=feval(f,x(i),y(i)); dy2=feval(f,x(i+1),y(i+1)); else y(i+1)=y(i)+h*(3*dy2-dy1)/2; P=feval(f,x(i+1),y(i+1)); y(i+1)=y(i)+h*(P+dy2)/2; dy1=dy2; dy2=feval(f,x(i+1),y(i+1)); end end A=[x,y]; 改进的4阶Adams预报校正系统： function A=CAdams4PC(f,a,b,N,ya) if N4 return; end h=(b-a)/N; x=zeros(1,N+1); y=zeros(1,N+1); x=a:h:b; y(1)=ya; for i=1:N if i4 k1=feval(f,x(i),y(i)); k2=feval(f,x(i)+h/2,y(i)+(h/2)*k1); k3=feval(f,x(i)+h/2,y(i)+(h/2)*k2); k4=feval(f,x(i)+h,y(i)+h*k3); y(i+1)=y(i)+(h/6)*(k1+2*k2+2*k3+k4); else if i==4 F=feval(f,x(i-3:i),y(i-3:i)); py=y(i)+(h/24)*(F*[-9,37,-59,55]); p=feval(f,x(i+1),py); F=[F(2)