此处参考教材为李荣华的《微分方程数值解法》
使用工具：Matlab
1. 算法
注: 最后一行应为k4,上面为笔误
2. 算法
I.需要求解的函数
function f=f1D(t,u,ft)
% ft为方程编号，u1D为精确解函数u(t),注意与f1D对应右端项函数f(t,u(t))
switch ft
case 1 %P10 习题1
%v=exp(-5*t);
f=-5*u;
case 2 % P27 例3.1
%v=(1+t.^2).^2;
f=4*t.*u.^(1/2);
case 3 %补充练习 P54 例1
%v=exp(-t);
f=-u;
end
II.对应精确解
function v=u1D(t,ft)
% ft为方程编号，u1D为精确解函数u（t）,注意与f1D 对应右端项函数f(t,u(t))
switch ft
case 1 %P10 习题1
v=exp(-5*t);
%f=-5*u;
case 2 % P27 例3.1
v=(1+t.^2).^2;
%f=4*t.*u.^(1/2);
case 3 %补充练习 P54 例1
%v=exp(-t);
f=-u;
end
III.求解问题的主函数
function ODESolver(ft,t0,T,h)
t=(t0:h:T);%t节点
n=length(t);%步数
ExactU=u1D(t,ft);%精确解
%--------------------------------------
% 1)向前欧拉，单步显式，u(i)=u(i-1)+h*f(t(i-1),u(i-1))
u=0.*t;%数值解初始化
u(1)=ExactU(1);%设初值
tic; %记录CPU计算开始时间
for i=2:n
u(i)=u(i-1)+h*f1D(t(i-1),u(i-1),ft);
end
time1=toc; %记录CPU计算结束时间
NumU1=u;%计算结束赋值
Err1=(NumU1(n)-ExactU(n));%T点处误差
RelErr1=Err1/ExactU(n);%相对误差
%--------------------------------------
% 2)向后Euler法，单步隐式，u(i)=u(i-1)+h*f(t(i),u(i))
u=0.*t;%数值解初始化
u(1)=ExactU(1);%设初值
tic; %记录CPU计算开始时间
for i=2:n
uc1=u(i-1)-1;%非线性迭代设初值
uc2=u(i-1);%非线性迭代设初值，因为迭代uc1=uc2,故初值uc2=u(i-1)
delta=10^(-5);%迭代阈值
while abs(uc2-uc1)>delta
uc1=uc2;
uc2=u(i-1)+h*f1D(t(i),uc1,ft);
end
u(i)=uc2;%非线性迭代结束
end
time2=toc; %记录CPU计算结束时间
NumU2=u; %计算结束赋值
Err2=NumU2(n)-ExactU(n); %T点处误差
RelErr2=Err2/ExactU(n); %相对误差
%--------------------------------------
% 3）梯形法，单步隐式法， u(i)=u(i-1)+(h/2)*(f(t(i-1),u(i-1))+f(t(i),u(i)))
u=0.*t; %数值解初始化
u(1)=ExactU(1); %设初值
tic; %记录CPU计算开始时间
for i=2:n
uc1=u(i-1)-1; %非线性迭代设置初值
uc2=u(i-1); %非线性迭代设置初值，因为迭代uc1=uc2,故初值uc2=u(i-1)
delta=10^(-5); %d迭代阈值
while abs(uc2-uc1)>delta
uc1=uc2;
uc2=u(i-1)+(h/2)*(f1D(t(i-1),u(i-1),ft)+f1D(t(i),uc1,ft));
end
u(i)=uc2; %非线性迭代结束
end
time3=toc; %记录CPU计算结束时间
NumU3=u; %计算结束赋值
Err3=NumU3(n)-ExactU(n); %T点处误差
RelErr3=Err3/ExactU(n); %相对误差
%--------------------------------------
% 4) 改进Euler法，单步显示法，uc（i）=u(i-1)+h*f(t(i-1),u(i-1)),u(i)=u(i-1)+(h/2)*(f(t(i-1),u(i-1))+f(t(i),uc(i)))
u=0.*t;%数值解法初始化
u(1)=ExactU(1);%设初值
tic;%记录CPU计算开始时间
for i=2:n
uc1=u(i-1)+h*4*t(i-1)*(u(i-1)^(1/2));%用向前Euler法预估
u(i)=u(i-1)+(h/2)*(f1D(t(i-1),u(i-1),ft)+f1D(t(i),uc1,ft));%用梯形法校正
end
time4=toc;%记录CPU计算结束时间
NumU4=real(u);%计算结束赋值
Err4=NumU4(n)-ExactU(n);%T点处误差
RelErr4=Err4/ExactU(n);%相对误差
%--------------------------------------
%5）二步法a=0,二步显式法,u(n+2)-(1+a)*u(n+1)+a*u(n)=(1/2)*h*((3-a)*f(n+1)-(1+a)*f(n))
u=0.*t;%数值解初始化
u(1)=ExactU(1);%设初值
u(2)=ExactU(2);%设初值
a=0;
tic;%记录CPU计算开始时间
for i=3:n
u(i)=(1+a)*u(i-1)-a*u(i-2)+(1/2)*h*((3-a)*f1D(t(i-1),u(i-1),ft)-(1+a)*f1D(t(i-2),u(i-2),ft));
end
time5=toc;%记录CPU计算结束时间
NumU5=u;%计算结束赋值
Err5=NumU5(n)-ExactU(n);%T点处误差
RelErr5=Err5/ExactU(n);%相对误差
%--------------------------------------
%6)二步法a=-5,二步显式法，u(n+2)-(1+a)*u(n+1)+a*u(n)=(1/2)*h*((3-a)*f(n+1)-(1+a)*f(n))
u=0.*t;%数值解初始化
u(1)=ExactU(1); %设初值
u(2)=ExactU(2); %设初值
a=-5;
tic; %记录CPU计算结束时间
for i=3:n
u(i)=(1+a)*u(i-1)-a*u(i-2)+(1/2)*h* ((3-a)*f1D(t(i-1),u(i-1),ft)-(1+a)*f1D(t(i-2),u(i-2),ft));
end
time6=toc; %记录CPU计算结束时间
NumU6=real(u); %计算结束赋值，由于u计算出现负数，开根号出现复数，故取实部
Err6=NumU6(n)-ExactU(n); %T点处误差
RelErr6=Err6/ExactU(n); %相对误差
%--------------------------------------
% 7）Simpson法，二步隐式法，u(i)=u(i-2)+(2*h/6)*(f(t(i-2),u(i-2))+4*f(t(i-1),u(i-1))+f(t(i),u(i)))
u=0.*t; %数值解初始化
u(1)=ExactU(1); %设初值
u(2)=ExactU(2); %设初值
tic; %记录CPU计算结束时间
for i=3:n
uc1=u(i-2); %非线性迭代设初值
uc2=u(i-1); %非线性迭代设初值，因为迭代uc1=uc2，故初值uc2=u(i-1)
delta=10^(-5); %迭代阈值，精度及运算时间受其影响
while abs(uc2-uc1)>delta
uc1=uc2;
uc2=u(i-2)+(h/3)*(f1D(t(i-2),u(i-2),ft)+4*f1D(t(i-1),u(i-1),ft)+f1D(t(i),uc1,ft));
end
u(i)=uc2; %非线性迭代结束
end
time7=toc; %记录CPU计算结束时间
NumU7=u; %计算结束赋值
Err7=NumU7(n)-ExactU(n);%T点处误差
RelErr7=Err7/ExactU(n);%相对误差
%--------------------------------------
% 8)四级四阶RK法，单步显式法
u=0.*t;%数值解初始化
u(1)=ExactU(1);%设初值
tic;%记录CPU开始计算的时间
for i=2:n
k1=f1D(t(i-1),u(i-1),ft);
k2=f1D(t(i-1)+h/2,u(i-1)+(h/2)*k1,ft);
k3=f1D(t(i-1)+h/2,u(i-1)+(h/2)*k2,ft);
k4=f1D(t(i-1)+h,u(i-1)+h*k3,ft);
u(i)=u(i-1)+(h/6)*(k1+2*k2+2*k3+k4);
end
time8=toc;%记录CPU计算结束时间
NumU8=u;%计算结束赋值
Err8=(NumU8(n)-ExactU(n));%T点处误差
RelErr8=Err8/ExactU(n);%相对误差
%--------------------------------------
%画图
figure;
plot(t,ExactU,'k-*',t,NumU1,'r-o',t,NumU2,'b-o',t,NumU3,'g-o',t,NumU4,'c-o',t,NumU5,'r-d',t,NumU7,'g-d',t,NumU8,'b-d');
legend('精确解','1)向前欧拉','2)向后欧拉','3)梯形法','4)改进Euler','5)二步法a=0','7)Simpson法','8）RK法','Location','NorthEastOutside');
%制表
f=figure;
data=([t;ExactU;NumU1;NumU2;NumU3;NumU4;NumU5;NumU6;NumU7;NumU8])';
data=[data;[0,0,Err1,Err2,Err3,Err4,Err5,Err6,Err7,Err8];[0,0,RelErr1,RelErr2,RelErr3,RelErr4,RelErr5,RelErr6,RelErr7,RelErr8];[0,0,time1,time2,time3,time4,time5,time6,time7,time8]];
colnames={'t','精确解','1)向前欧拉','2)向后欧拉','3)梯形法','4)改进Euler','5)二步法a=0','6)二步法a=-5','7)Simpson法','8）RK法'};
tab=uitable(f,'Data',data,'ColumnName',colnames,'Position',[50 50 800 500]);
3.给出简单结果
[^此处图略，仅给出简单数字结果；感兴趣的读者可以复制代码运行一下]
P10页：习题1
（1） 步长0.1：
在命令行窗口输入
>>ODESolver(1,0,1,0.1)
所得最后一步误差结果：
Err1 =-0.0058
Err2 =0.0106
Err3 =-6.9002e-04
Err4 =-0.0098
Err5 =0.0045
Err6 =-2.1480e+04
Err7 =-3.0847e-04
Err8 =2.6728e-05
（2）步长0.05：
在命令行窗口输入
>>ODESolver(1,0,1,0.05)
所得结果：
Err1 =-0.0036
Err2 =0.0048
Err3 =-1.7339e-04
Err4 =-0.0069
Err5 = 9.8516e-04
Err6 = -2.1437e+10
Err7 =-1.2741e-05
Err8 = 1.3516e-06
P27页：例题3.1
（1）步长0.1：
在命令行窗口输入
>>ODESolver(2,0,5,0.1)
所得结果：
Err1 = -3.5381
Err2 = 4.0508
Err3 = 0.0806
Err4 = -0.0763
Err5 = -0.3675
Err6 = -6.9711e+08
Err7 = -1.3874e-06
Err8 = -9.5303e-05
（2）步长0.05：
在命令行窗口输入
>>ODESolver(2,0,2,0.05)
所得结果：
Err1 = -1.8298
Err2 = 1.9580
Err3 = 0.0201
Err4 = -0.0196
Err5 = -0.0963
Err6 = -5.4646e+21
Err7 = -1.4567e-06
Err8 = -6.2349e-06
P54 例1
（1）步长0.1：
在命令行窗口输入
>>ODESolver(3,0,2,0.1)
所得结果：
Err1 = -0.0138
Err2 = 0.0133
Err3 = -2.2421e-04
Err4 = -0.0974
Err5 = 0.0011
Err6 = -1.2434e+08
Err7 = 1.8151e-08
Err8 = 2.4519e-07
4.对算例结果进行简单的分析
从整体来看：
第一道题在靠近1的附近，数值解精确度都比较高。第二道题在靠近0的附近，数值解精确度都比较高。其中，向前Euler法和向后Euler法落在精确解的两侧。图表的最后三行分别表示最后一步的数值误差、最后一步的相对误差、每个算法所花费的时间。
从步长来看：
图像中可以较直观的发现，h=0.05比h=0.1的精度整体更高一些。当步长变小时，对结果的分析更加清晰。
从八种不同的算法来看：
利用直观的图像以及数据表，Simpson方法误差最小，其次是RK法，二步法的误差最大（因为这种算法不稳定，可从图表中看出当a=-5时，二步法显示出不稳定性，且当步长缩小，时间增加时，不稳定性表现得非常明显）。在所有的单步法中，梯形法和改进欧拉法效果也比较好。
1282
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
