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

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

d x d t = f ( x , t ) , x ( t 0 ) = x 0 \frac{d x}{d t}=f(x, t), \quad x\left(t_{0}\right)=x_{0} dtdx​=f(x,t),x(t0​)=x0​

1. 前向欧拉法

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

x ( t 0 + h ) = x ( t 0 ) + h x ′ ( t 0 ) + 1 2 x ′ ′ ( ξ ) h 2 , t 0 < ξ < t 0 + h 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}, \quad t_{0}<\xi<t_{0}+h x(t0​+h)=x(t0​)+hx′(t0​)+21​x′′(ξ)h2,t0​<ξ<t0​+h

x k + 1 = x k + h x k ′ + O ( h 2 ) = x k + h f ( x k , t k ) + O ( h 2 ) x_{k+1}=x_{k}+h x'_k+O\left(h^{2}\right)=x_{k}+h f\left(x_{k}, t_{k}\right)+O\left(h^{2}\right) xk+1​=xk​+hxk′​+O(h2)=xk​+hf(xk​,tk​)+O(h2)

2. 改进的欧拉法

在原来前向欧拉法的基础上泰勒展开使用了前面两项：

x k + 1 = x k + h x k ′ + 1 2 h 2 x k ′ ′ + O ( h 3 ) x_{k+1}=x_{k}+h x^{\prime}_k+\frac{1}{2} h^{2} x_{k}^{\prime \prime}+O\left(h^{3}\right) xk+1​=xk​+hxk′​+21​h2xk′′​+O(h3)

这里使用：

x k ′ ′ = x k + 1 ′ − x k ′ h x_{k}^{\prime \prime}=\frac{x_{k+1}^{\prime}-x_{k}^{\prime}}{h} xk′′​=hxk+1′​−xk′​​

于是我们有：

x k + 1 = x k + h 2 ( x k ′ + x k + 1 ′ ) + O ( h 3 ) x_{k+1}=x_{k}+\frac{h}{2}\left(x_{k}^{\prime}+x_{k+1}^{\prime}\right)+O\left(h^{3}\right) xk+1​=xk​+2h​(xk′​+xk+1′​)+O(h3)

也就是：

x k + 1 = x k + h 2 [ f ( x k , t k ) + f ( x k + 1 , t k + 1 ) ] + O ( h 3 ) 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) xk+1​=xk​+2h​[f(xk​,tk​)+f(xk+1​,tk+1​)]+O(h3)

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

对比前向欧拉法，改进欧拉法的右边不使用 x k + 1 x_{k+1} xk+1​（我们还不知道），但是我们可以用前向欧拉法计算的 x k + h f ( x k , t k ) x_{k}+h f\left(x_{k}, t_{k}\right) xk​+hf(xk​,tk​)来代替 x k + 1 x_{k+1} xk+1​，于是我们有

x k + 1 = x k + 1 2 ( k 1 + k 2 ) + O ( h 3 ) , 其中： k 1 = h f ( x k , t k ) , k 2 = h f ( x k + h , t k + k 1 ) x_{k+1}=x_{k}+\frac{1}{2}\left(k_{1}+k_{2}\right)+O\left(h^{3}\right), \\\text{其中：}k_{1}=h f\left(x_{k}, t_{k}\right), k_{2}=h f\left(x_{k}+h, t_{k}+k_{1}\right) xk+1​=xk​+21​(k1​+k2​)+O(h3),其中：k1​=hf(xk​,tk​),k2​=hf(xk​+h,tk​+k1​)

对比一下前向欧拉法：

x k + 1 = x k + k 1 + O ( h 2 ) , k 1 = h f ( x k , t k ) x_{k+1}=x_{k}+k_{1}+O\left(h^{2}\right), \quad k_{1}=h f\left(x_{k},t_{k} \right) xk+1​=xk​+k1​+O(h2),k1​=hf(xk​,tk​)

例子

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

clcclear% 测试三个不同的步长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;% 设置初始条件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);% 定义微分方程f = inline('xx+tt','tt','xx');for k = 1:nt(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));enddiff1_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_timessubplot(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('t','Fontsize',18)ylabel('x','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_timessemilogy(tplot_res{i},diff1_res{i},'b-',tplot_res{i},diff2_res{i},'r:')hold onnum1 = 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')endxlabel('t','Fontsize',18)ylabel('|relative error|','Fontsize',18)legend({'Euler method','modified Euler method'},'Location','best')title('Two Euler method''s relative error')legend boxoff;

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