matlab解常微分方程常用数值解法1：前向欧拉法和改进的欧拉法

总结和记录一下matlab求解常微分方程常用的数值解法，本文先从欧拉法和改进的欧拉法讲起。

$$\frac{dx}{dt}=f(x,t),x\left(t_0\right)=x_0$$

1. 前向欧拉法

前向欧拉法使用了泰勒展开的第一项线性项逼近。

$$x\left(t_0+h\right)=x\left(t_0\right)+h{x}^{\prime }\left(t_0\right)+\frac{1}{2}{x}^{\prime \prime }(\xi )h^2,t_0<\xi <t_0+h$$

$$x_{k+1}=x_k+h{x}_k^{\prime }+O\left(h^2\right)=x_k+hf\left(x_k,t_k\right)+O\left(h^2\right)$$

2. 改进的欧拉法

在原来前向欧拉法的基础上泰勒展开使用了前面两项：
$$x_{k+1}=x_k+h{x}_k^{\prime }+\frac{1}{2}{h}^2{x}_k^{\prime \prime }+O\left(h^3\right)$$

这里使用：

$$x_k^{\prime \prime }=\frac{x_{k+1}^{\prime }-x_k^{\prime }}{h}$$

于是我们有：

$$x_{k+1}=x_k+\frac{h}{2}\left(x_k^{\prime }+x_{k+1}^{\prime }\right)+O\left(h^3\right)$$

也就是：

$$x_{k+1}=x_k+\frac{h}{2}\left[f\left(x_k,t_k\right)+f\left(x_{k+1},t_{k+1}\right)\right]+O\left(h^3\right)$$

我们怎么计算$f(x_{k+1},t_{k+1})$呢，因为我们还不知道$x_{k+1}$。

对比前向欧拉法，改进欧拉法的右边不使用$x_{k+1}$（我们还不知道），但是我们可以用前向欧拉法计算的$x_k+hf\left(x_k,t_k\right)$来代替$x_{k+1}$，于是我们有
$$\begin{gathered} x_{k+1}=x_k+\frac{1}{2}\left(k_1+k_2\right)+O\left(h^3\right), \\ \text{其中：}{k}_1=hf\left(x_k,t_k\right),k_2=hf\left(x_k+h,t_k+k_1\right) \end{gathered}$$

对比一下前向欧拉法：

$$x_{k+1}=x_k+k_1+O\left(h^2\right),k_1=hf\left(x_k,t_k\right)$$

例子

$$x^{\prime }=x+t,x(0)=1$$

clear% 测试三个不同的步长
test_times = 3;
% 保存时间、差分时间和步长
h_res=ones(1,test_times);
t_res=cell(1,test_times);%时间
tplot_res=cell(1,test_times);%差分的时间，比时间长度少1
% 保存两种数值方法和解析解的计算结果
x_euler_res=cell(1,test_times);
x_modified_res=cell(1,test_times);
x_exact_res=cell(1,test_times);
% 保存误差
diff1_res=cell(1,test_times);
diff2_res=cell(1,test_times);for i = 1:test_times
% 设置步长间隔和步长数
h = 1/10^i; n = 10/h;
% set up initial conditions
t=zeros(n+1,1); t(1) = 0;
x_euler=zeros(n+1,1); x_euler(1) = 1;
x_modified=zeros(n+1,1); x_modified(1) = 1;
x_exact=zeros(n+1,1); x_exact(1) = 1;
% 设置不同的比较误差的图
diff1 = zeros(n,1); diff2 = zeros(n,1); tplot = zeros(n,1);
% define right side of differential equation, Equation 1.7.10
f = inline('xx+tt','tt','xx');
for k = 1:n
t(k+1) = t(k) + h;
% 计算解析解
x_exact(k+1) = 2*exp(t(k+1)) - t(k+1) - 1;
% 使用前向欧拉法计算
k1 = h * f(t(k),x_euler(k));
x_euler(k+1) = x_euler(k) + k1;
tplot(k) = t(k+1);
diff1(k) = x_euler(k+1) - x_exact(k+1);
diff1(k) = abs(diff1(k) / x_exact(k+1));
% 使用改进欧拉法计算
k1 = h * f(t(k),x_modified(k));
k2 = h * f(t(k+1),x_modified(k)+k1);
x_modified(k+1) = x_modified(k) + 0.5*(k1+k2);
diff2(k) = x_modified(k+1) - x_exact(k+1);
diff2(k) = abs(diff2(k) / x_exact(k+1));
end
diff1_res{i}=diff1;
diff2_res{i}=diff2;
tplot_res{i}=tplot;
h_res(i)=h;
x_euler_res{i}=x_euler;
x_modified_res{i}=x_modified;
x_exact_res{i}=x_exact;
t_res{i}=t;
endfigure
for i=1:test_times
subplot(2,2,i)
plot(t_res{i},x_exact_res{i},'k-',t_res{i},x_euler_res{i},'b-',t_res{i},x_modified_res{i},'r:')
xlabel('TIME','Fontsize',18)
ylabel('|RELATIVE ERROR|','Fontsize',18)
legend({'Analytical method','Euler method','modified Euler method'},'Location','best')
legend boxoff;
title(['h = ',num2str(h_res(i))]);
end
subplot(2,2,4)% 计算相对误差
for i=1:test_times
semilogy(tplot_res{i},diff1_res{i},'b-',tplot_res{i},diff2_res{i},'r:')
hold on
num1 = 0.2*10/h_res(i); num2 = 0.8*10/h_res(i);
text(3,diff1_res{i}(num1),['h = ',num2str(h_res(i))],'Fontsize',12,...
'HorizontalAlignment','right',...
'VerticalAlignment','bottom')
text(9,diff2_res{i}(num2),['h = ',num2str(h_res(i))],'Fontsize',12,...
'HorizontalAlignment','right',...
'VerticalAlignment','bottom')
end
xlabel('TIME','Fontsize',18)
ylabel('|RELATIVE ERROR|','Fontsize',18)
legend({'Euler method','modified Euler method'},'Location','best')
legend boxoff;

我们对各个不同的步长进行了比较，并比较了它们的误差，发现改进的欧拉法要比前向欧拉法的精度更高。随着步长的变小，误差也在变小。