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

本文将使用 Matlab 编程实现隐式 Euler 法，求解如下常微分方程初值问题：

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

首先，我们需要将常微分方程转化为离散形式，隐式 Euler 法的离散形式为：

u(i+1) = u(i) + h*(-4*u(i+1) + 6*v(i+1))
v(i+1) = v(i) + h*(3*u(i+1) - 7*v(i+1))

其中，h 为时间步长，i 为时间步数。将初始值代入，我们可以得到：

u(1) = 1
v(1) = 2

接下来，我们可以用循环来求解：

h = 0.01; % 时间步长
t = 0:h:1; % 时间步数
N = length(t); % 时间步数总数

u = zeros(N, 1); % 初始化 u 数组
v = zeros(N, 1); % 初始化 v 数组

u(1) = 1; % 初始值
v(1) = 2; % 初始值

for i = 2:N
    % 隐式 Euler 法
    u(i) = (u(i-1) + h*6*v(i-1)) / (1 + 4*h);
    v(i) = (v(i-1) + h*3*u(i)) / (1 + 7*h);
end

接下来，我们可以计算真实解，并绘制误差图：

u_true = 2*exp(-2*t) + 3*exp(-t); % 真实解
v_true = -exp(-2*t) - exp(-t); % 真实解

u_error = abs(u - u_true'); % u 数值解误差
v_error = abs(v - v_true'); % v 数值解误差

plot(t, u_error, 'r', t, v_error, 'b'); % 绘制误差图
xlabel('t');
ylabel('Numerical Solution Error');
legend('u', 'v');

最终的图像如下所示：

可以看到，隐式 Euler 法的误差比较小，说明该方法比较有效。