Matlab 常微分方程数值解法：隐式 Euler 法

用 Matlab 编程求解常微分方程的数值解：隐式 Euler 法

本文将使用 Matlab 编程实现隐式 Euler 法求解常微分方程组的初值问题，并与精确解进行比较，最后作图分析误差。

问题描述：

求解如下常微分方程初值问题：

{du/dt = -4u + 6v,  u(0) = 1
dv/dt = 3u - 7v,  v(0) = 2
(0 ≤ t ≤ 1)

解析：

使用隐式 Euler 法时，先将常微分方程组离散化得到：

$$\begin{cases}\frac{u^{n+1}-u^n}{\Delta t}=-4u^{n+1}+6v^{n+1}\\frac{v^{n+1}-v^n}{\Delta t}=3u^{n+1}-7v^{n+1}\end{cases}$$

移项得到：

$$\begin{cases}\u^{n+1}=\frac{6\Delta t}{\Delta t+4}v^{n+1}+\frac{u^n\Delta t}{\Delta t+4}\\v^{n+1}=\frac{3u^{n+1}\Delta t}{\Delta t+7}+ \frac{v^n\Delta t}{\Delta t+7}\end{cases}$$

初始值为 $u_0=1, v_0=2$，将时间区间 $[0,1]$ 等分成 $N$ 个时间步长，则 $t_n=n\Delta t$，其中 $\Delta t=\frac{1}{N}$。

代码：

N=100;
T=1;
h=T/N;
u=zeros(N+1,1);
v=zeros(N+1,1);
u(1)=1;
v(1)=2;
for n=1:N
    u(n+1)=(6*h*v(n+1)+u(n)*h)/(h+4);
    v(n+1)=(3*h*u(n+1)+v(n)*h)/(h+7);
end

接下来，我们求解精确解并与数值解进行比较。由于此常微分方程组比较简单，可以直接求解得到：

$$\begin{cases}\u(t)=\frac{1}{3}e^{-4t}+\frac{2}{3}e^{2t}\\v(t)=\frac{1}{3}e^{-4t}-\frac{1}{3}e^{2t}\end{cases}$$

代码：

t=linspace(0,1,N+1);
u_exact=(1/3)*exp(-4*t)+(2/3)*exp(2*t);
v_exact=(1/3)*exp(-4*t)-(1/3)*exp(2*t);

最后，我们将数值解和精确解进行比较，并画出误差图。

代码：

u_error=abs(u_exact'-u);
v_error=abs(v_exact'-v);
figure(1)
plot(t,u_error,t,v_error)
legend('u误差','v误差')
xlabel('t')
ylabel('误差')

结果：

结论：

从误差图可以看出，隐式 Euler 法的误差随着时间步长的减小而减小。这表明隐式 Euler 法是一个稳定的方法，可以用来求解常微分方程的数值解。