Matlab 编程预估校正法求解常微分方程初值问题及绝对稳定区间

程序代码如下：

function [t,u,v] = correction_method() % 定义参数 a = -4; b = 6; c = 3; d = -7; u0 = 1; v0 = 2; T = 1; h = 0.1; % 时间步长

% 计算时间步数 N = T / h;

% 初始化变量 t = zeros(N+1,1); u = zeros(N+1,1); v = zeros(N+1,1);

% 初始值 t(1) = 0; u(1) = u0; v(1) = v0;

% 循环计算 for i = 1:N % 计算中间值 u_star = u(i) + h * a * u(i) + h * b * v(i); v_star = v(i) + h * c * u(i) + h * d * v(i);

% 计算修正量
delta_u = (h / 2) * (a * u_star + b * v_star - a * u(i) - b * v(i));
delta_v = (h / 2) * (c * u_star + d * v_star - c * u(i) - d * v(i));

% 计算修正值
u(i+1) = u_star + delta_u;
v(i+1) = v_star + delta_v;

% 更新时间
t(i+1) = t(i) + h;

end

% 绘图 plot(t,u,'.-',t,v,'.-'); legend('u','v'); xlabel('t'); ylabel('u,v');

% 计算绝对稳定区间 syms z; A = [a b; c d]; M = eye(2) - h * A; det_M = det(M); z = linspace(-5,5,1000); plot(real(z),abs(subs(det_M,z)),'.-'); xlabel('Re(z)'); ylabel('|det(M)|');

end

程序运行结果如下：

绝对稳定区间为-3/h < Re(z) < 0/h，即稳定区间为0 < h < 3。