常微分方程初值问题数值解法报告.docVIP

常微分方程初值问题数值解法 081180081 唐方爽 电子 题：选用一种数值解法，求解初值问题取h=0.1，0.01，0.001，……，计算x=1处解的近似值。可用方法：（1）梯形法；（2）中点法或修正中点法；（3）标准R-K法；（4）修正哈明预测-校正法。 解： 梯形法： h=0.1; x=0:h:1; y=zeros(size(x)); y(1)=1; for i = 1:length(x)-1 y(i+1)=y(i)+h*f(x(i),y(i)); y(i+1)=y(i)+h/2*(f(x(i),y(i))+f(x(i+1),y(i+1))); end figure,semilogy(x,y,o),hold on,title(梯形法) 实验结果如下： h=0.1时，x=1处的近似解为26.3841 h=0.01时，x=1处的近似解为218.7920 h=0.001时，x=1处的近似解为2.0993e+003 中点法： clc,close all h=0.1; x=0:h:1; y=zeros(size(x)); y(1)=1; y(2)=y(1)+h*f(x(1),y(1)); for i = 2:length(x)-1 y(i+1)=y(i-1)+2*h*f(x(i),y(i)); end y(end) figure,semilogy(x,y,o),hold on,title(中点法) h=0.1时，x=1处的近似解为15.1256 h=0.01时，x=1处的近似解为134.2798 h=0.001时，x=1处的近似解为1.3077e+003 标准R-K法： clc,close all h=0.1; x=0:h:1; y=zeros(size(x)); y(1)=1; for i = 1:length(x)-1 K1=h*f(x(i),y(i)); K2=h*f(x(i)+h/2,y(i)+K1/2); K3=h*f(x(i)+h/2,y(i)+K2/2); K4=h*f(x(i)+h,y(i)+K3); y(i+1)=y(i)+(K1+2*K2+2*K3+K4)/6; end y(end) figure,semilogy(x,y,o),hold on,title(标准R-K法) h=0.1时，x=1处的近似解为85.3348 h=0.01时，x=1处的近似解为823.4736 h=0.001时，x=1处的近似解为8.2027e+003 修正哈明预测校正方法： h=0.1; x=0:h:1; y=zeros(size(x)); p=zeros(size(x)); c=zeros(size(x)); m=zeros(size(x)); y(1)=1; for i = 1:3 %用标准龙格库塔法算出前三点 K1=h*f(x(i),y(i)); K2=h*f(x(i)+h/2,y(i)+K1/2); K3=h*f(x(i)+h/2,y(i)+K2/2); K4=h*f(x(i)+h,y(i)+K3); y(i+1)=y(i)+(K1+2*K2+2*K3+K4)/6; end for i = 4:length(x)-1 p(i+1)=y(i-3)+4/3*h*(2*f(x(i),y(i))-f(x(i-1),y(i-1))+2*f(x(i-2),y(i-2))); m(i+1)=p(i+1)+112/121*(c(i)-p(i)); c(i+1)=(9*y(i)-y(i-2))/8+3/8*h*(f(x(i+1),m(i+1))+2*f(x(i),y(i))-f(x(i-1),y(i-1))); y(i+1)=c(i+1)-9/121*(c(i+1)-p(i+1)); end y(end) figure,semilogy(x,y,o),hold on,title(修正哈明预测-校正法) h=0.1时，x=1处的近似解为46.2812 h=0.01时，x=1处的近似解为444.9971 h=0.001时，x=1处的近似解为4.4312e+003 X=1处的精确解是无穷大，因此x=1处数值解越大越精确 可见，各种方法中，相同步长条件下，标准R-K法精度最高，其次是修正哈明预测校正法，再次是梯形法，最差是中点法