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

1 常微分方程的差分方法实验实验目的及要求：（1）深入理解常微分方程的差分方法的原理，学会用差分方法解决某些实际的常微分方程问题，比较这些方法解题的不同之处2）熟悉 Matlab 编程环境，利用 Matlab 实现具体的常微分方程用 Matlab 软件实现欧拉方法、 改进的欧拉方法、 龙格-库塔方法和亚当姆斯方法，并用实例在计算机上计算实验内容 ：实验题目 （1）取 h=0.1，用欧拉方法、改进的欧拉方法、四阶龙格-库塔方法求解初值问题：0)0(2112 2yyxy40x并与精确解221xxy比较计算结果 . 2 （2）分别用二阶亚当姆斯预估校正系统、改进的四阶亚当姆斯预估校正系统求解 初值问题0)0(1yyy，，取181.0)2.0(2.0yh，，计算)0.1(y实验过程 ：(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); 1. Euler 方法：首先编写.Euler m文件，如下：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']; 1. 改进的 Euler 方法：编写文件MendEulerm文件，如下：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 3 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']; 1. 四阶 Runge-Kutta 方法：编写4.Rungkutam文件，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; 4 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 N> f=@f2;a=0;b=4;N=40;ya=0; >> E=Euler(f,a,b,N,ya); >> y=solvef2(a:(b-a)/N:b); >> m=[E,y'] m = 0 0 0 0.10000000000000 0.10000000000000 0.09900990099010 0.20000000000000 0.19700990099010 0.19230769230769 0.30000000000000 0.28540116692632 0.27522935779817 0.40000000000000 0.36085352097579 0.34482758620690 0.50000000000000 0.42101736480739 0.40000000000000 0.60000000000000 0.46556624051352 0.44117647058824 0.70000000000000 0.49574526741705 0.46979865771812 6 0.80000000000000 0.51370668734350 0.48780487804878 0.90000000000000 0.52190338497531 0.49723756906077 1.00000000000000 0.52267537511010 0.50000000000000 1.10000000000000 0.51803746556081 0.49773755656109 1.20000000000000 0.50961377119415 0.49180327868852 1.30000000000000 0.49865613859339 0.48327137546468 1.40000000000000 0.48609927087160 0.47297297297297 1.50000000000000 0.47262455442701 0.46153846153846 1.60000000000000 0.45871899130677 0.44943820224719 1.70000000000000 0.44472425635012 0.43701799485861 1.80000000000000 0.43087526438692 0.42452830188679 1.90000000000000 0.41732947135520 0.41214750542299 2.00000000000000 0.40418866779251 0.40000000000000 2.10000000000000 0.39151497195813 0.38817005545287 2.20000000000000 0.37934246565956 0.37671232876712 2.30000000000000 0.36768561208025 0.36565977742448 2.40000000000000 0.35654532140646 0.35502958579882 2.50000000000000 0.34591330757137 0.34482758620690 2.60000000000000 0.33577520774866 0.33505154639175 2.70000000000000 0.32611280765907 0.32569360675513 2.80000000000000 0.31690562117133 0.31674208144796 2.90000000000000 0.30813200381990 0.30818278427205 3.00000000000000 0.29976993002539 0.30000000000000 3.10000000000000 0.29179752783591 0.29217719132893 3.20000000000000 0.28419343907371 0.28469750889680 3.30000000000000 0.27693705406423 0.27754415475189 3.40000000000000 0.27000865661335 0.27070063694268 3.50000000000000 0.26338950512361 0.26415094339623 3.60000000000000 0.25706186865308 0.25787965616046 3.70000000000000 0.25100903157223 0.25187202178353 3.80000000000000 0.24521527672616 0.24611398963731 3.90000000000000 0.23966585427601 0.24059222702036 4.00000000000000 0.23434694139690 0.23529411764706 >> 2. 改进的 Euler 方法： >> f=@f2;a=0;b=4;N=40;ya=0; >> E=MendEuler(f,a,b,N,ya); >> y=solvef2(a:(b-a)/N:b); >> m=[E,y'] m = 0 0 0 7 0.10000000000000 0.09850495049505 0.09900990099010 0.20000000000000 0.19129157451772 0.19230769230769 0.30000000000000 0.27373370103545 0.27522935779817 0.40000000000000 0.34293131547441 0.34482758620690 0.50000000000000 0.39782199226245 0.40000000000000 0.60000000000000 0.43885373990985 0.44117647058824 0.70000000000000 0.46746146344897 0.46979865771812 0.80000000000000 0.48555880905455 0.48780487804878 0.90000000000000 0.49515605455715 0.49723756906077 1.00000000000000 0.49812534779304 0.50000000000000 1.10000000000000 0.49608671330636 0.49773755656109 1.20000000000000 0.49037501355228 0.49180327868852 1.30000000000000 0.48205289737796 0.48327137546468 1.40000000000000 0.47194514636326 0.47297297297297 1.50000000000000 0.46067950299100 0.46153846153846 1.60000000000000 。