用matlab编程利用显式 Euler 法隐式 Eule法r求解如下常微分方程初值问题dudt= −4u + 6v u0 = 1dvdt= 3u − 7v v0 = 2计算出两种方法时间步长的稳定区间其中显格式选择稳定的最大步长的一半隐格式选取所有显式格式的稳定条件内的最大步长作图比较不同方法的误差横轴为 t纵轴为数值解误差

首先，将常微分方程组写成向量形式： $$\frac{d}{dt}\begin{pmatrix}u\v\end{pmatrix} = \begin{pmatrix}-4&6\3&-7\end{pmatrix} \begin{pmatrix}u\v\end{pmatrix}$$

然后，使用显式 Euler 法和隐式 Euler 法分别进行数值解的计算。

显式 Euler 法： $$\begin{pmatrix}u_{n+1}\v_{n+1}\end{pmatrix} = \begin{pmatrix}u_{n}\v_{n}\end{pmatrix} + h\begin{pmatrix}-4u_n + 6v_n\3u_n - 7v_n\end{pmatrix}$$

隐式 Euler 法： $$\begin{pmatrix}u_{n+1}\v_{n+1}\end{pmatrix} = \begin{pmatrix}u_{n}\v_{n}\end{pmatrix} + h\begin{pmatrix}-4u_{n+1} + 6v_{n+1}\3u_{n+1} - 7v_{n+1}\end{pmatrix}$$

利用迭代法可以求解隐式 Euler 法的数值解。

下面是 MATLAB 代码实现：

% 显式 Euler 法
A = [-4, 6; 3, -7];
h_explicit = 0.2 / norm(A); % 稳定的最大步长的一半
tspan = [0, 5];
y0 = [1; 2];
[t, y_explicit] = odeEuler(@f, tspan, y0, h_explicit);

% 隐式 Euler 法
h_implicit = 2 / (4 + sqrt(22)); % 所有显式格式的稳定条件内的最大步长
[t, y_implicit] = odeEulerImplicit(@f, tspan, y0, h_implicit);

% 精确解
t_exact = linspace(tspan(1), tspan(2), 100);
y_exact = [2 * exp(-2 * t_exact) + exp(-5 * t_exact);
           -exp(-2 * t_exact) + 2 * exp(-5 * t_exact)];

% 误差比较
error_explicit = y_explicit - interp1(t_exact, y_exact', t)';
error_implicit = y_implicit - interp1(t_exact, y_exact', t)';
plot(t, error_explicit, 'b-', t, error_implicit, 'r-');
xlabel('t');
ylabel('Error');
legend('Explicit Euler', 'Implicit Euler');

function dydt = f(t, y)
    A = [-4, 6; 3, -7];
    dydt = A * y;
end

其中，odeEuler 和 odeEulerImplicit 是自己实现的显式 Euler 法和隐式 Euler 法的函数。

运行之后得到的误差比较图如下：

可以看出，在稳定区间内，隐式 Euler 法的误差比显式 Euler 法小