微分方程数值解法（1）——常微分方程初值问题的数值解法

此处参考教材为李荣华的《微分方程数值解法》
使用工具：Matlab

1. 算法

注: 最后一行应为k₄,上面为笔误

2. 算法

I.需要求解的函数

function f=f1D(t,u,ft)% ft为方程编号，u1D为精确解函数u(t),注意与f1D对应右端项函数f(t,u(t))switch ftcase 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 ftcase 1 %P10 习题1v=exp(-5*t);%f=-5*u;case 2 % P27 例3.1v=(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:nu(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:nuc1=u(i-1)-1;%非线性迭代设初值uc2=u(i-1);%非线性迭代设初值，因为迭代uc1=uc2,故初值uc2=u(i-1)delta=10^(-5);%迭代阈值while abs(uc2-uc1)>deltauc1=uc2;uc2=u(i-1)+h*f1D(t(i),uc1,ft);endu(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:nuc1=u(i-1)-1; %非线性迭代设置初值uc2=u(i-1); %非线性迭代设置初值，因为迭代uc1=uc2,故初值uc2=u(i-1)delta=10^(-5); %d迭代阈值while abs(uc2-uc1)>deltauc1=uc2;uc2=u(i-1)+(h/2)*(f1D(t(i-1),u(i-1),ft)+f1D(t(i),uc1,ft));endu(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:nuc1=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:nuc1=u(i-2); %非线性迭代设初值uc2=u(i-1); %非线性迭代设初值，因为迭代uc1=uc2，故初值uc2=u(i-1)delta=10^(-5); %迭代阈值，精度及运算时间受其影响while abs(uc2-uc1)>deltauc1=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));endu(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:nk1=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]);