Matlab龙格-库塔法作业：P、Z求解及误差绘图故障排查

我在Matlab课程作业中有两项任务：
a) 求解微分方程组：
$$\frac{dP}{dt} = rP(1 - \frac{P^2}{\alpha2+P^2}) - RmZ$$
$$\frac{dZ}{dt} = \gamma RmZ \frac{P}{K+P} - \mu Z$$
得到P与Z的数值解；
b) 通过计算(|P(h/2)-P(h)|)（h为步长）绘制P和Z的误差函数。

目前任务a的结果正确，但任务b无法得到正确结果，请求排查代码问题。

相关代码 #

主程序 ah_alg.m

% File for main program (ah_alg)
clear all
close all 
clc

h = 0.5;                  % Step length
h2 = h/2
tspan = 0:h:400;       %Start value, time expressed as days

% Initial conditions
p0 = 20
z0 = 5
x0 = [p0,z0]'; % Write everything in vector form
r = 0.3
K = 108
R = 0.7


%Parameters
alpha = 5.7
mu = 0.024
gamma = 0.05

% Solving Runge knutta 4

X(:,1) = x0;
xin = x0;
xin2 = x0;
Err = [0; 0]; % Matrix for saving error for p and x

time = tspan(1);

for i=1:tspan(end)/h
    
    pout = rk4(@(t,x)lorenz(t,x,K,R,r,alpha,mu,gamma),h,time,xin);
    pout2 = rk4(@(t,x)lorenz(t,x,K,R,r,alpha,mu,gamma),h2,time,xin2);
    X = [X pout];
    if time <= 200
      Err(1,i) = abs(pout2(1)-pout(1));  %Error for p
      Err(2,i) = abs(pout2(2)-pout(2)); %Error for z
    end
    xin = pout;
    xin2 = pout2;
    time = time+h;
end
pout
Err(:,(201:400))= [];
length(Err(1,:))
%length(tspan)
plot(log(tspan(1,(1:201))),log(Err(1,:))) % Erf for p
xlabel('logtime');
%hold on
%plot(log(tspan(1,:)),log(Err(2,:)))
%hold on
%plot(tspan, X(1,:));
%hold on
%plot(tspan,X(2,:))

微分方程组函数 lorenz.m

% Lorenz function
function dp = lorenz(t,x,K,R,r,alpha,mu,gamma)
  dp = [
  r*x(1)*(1-x(1)/K) - R*x(2)*(x(1)^2/(alpha^2 + x(1)^2));
  gamma*R*x(2)*(x(1)^2/(alpha^2 + x(1)^2)) - mu*x(2);
  ];

RK4求解函数 rk4.m

% Rk4
function pout = rk4(f,h,tk,xk)
%Rungeknutta 4 method
%   This function calculates Rk4
% Startvalue for k1, t = t0

  k1 = f(tk,xk);  
  k2 = f(tk+h/2,xk+h*k1/2);
  k3 = f(tk+h/2,xk+h*k2/2);
  k4 = f(tk+h,xk+h*k3);
  pout = xk + (h/6)*(k1+k2+k3+k4);

问题分析与修复方案 #

1. 核心错误：步长迭代的时间不匹配

任务b要求对比同一时间点下h步长和h/2步长的解，但原代码中：

• 每次循环，xin走1步h，时间推进h；
• xin2只走1步h/2，时间仅推进h/2；
  导致两者的时间始终相差h/2，对比的不是同一时刻的P/Z值，误差计算完全错误。

2. 次要错误：微分方程组实现不符题目要求

原lorenz函数的方程和题目给出的不一致：

• 题目中dP/dt的第二项是(-RmZ)，代码中多了(P^2/(\alpha2+P^2))因子；
• 题目中dZ/dt的第一项是(\gamma RmZ \frac{P}{K+P})，代码中把分母写成(\alpha^2+P2)，分子写成(P^2)。

3. 其他问题：绘图与矩阵初始化的细节错误

• Err初始化为[0;0]，循环中动态扩展容易出现索引混乱；
• 绘图时使用tspan(1,(1:201))包含t=0，log(0)会产生数值错误。

修复后的代码 #

修正后的lorenz.m

% Lorenz function（匹配题目方程组）
function dp = lorenz(t,x,K,R,r,alpha,mu,gamma)
    dp = [
        r*x(1)*(1 - x(1)^2/(alpha^2 + x(1)^2)) - R*x(2);
        gamma*R*x(2)*(x(1)/(K + x(1))) - mu*x(2);
    ];
end

修正后的主程序ah_alg.m

% File for main program (ah_alg)
clear all
close all 
clc

h = 0.5;                  
h2 = h/2;
tspan = 0:h:400;          
tspan_h2 = 0:h2:400;      % h/2步长的完整时间轴

% 初始条件与参数
p0 = 20;
z0 = 5;
x0 = [p0,z0]'; 
r = 0.3;
K = 108;
R = 0.7;
alpha = 5.7;
mu = 0.024;
gamma = 0.05;

% 预分配矩阵存储解
X_h = zeros(2, length(tspan));
X_h(:,1) = x0;
X_h2 = zeros(2, length(tspan_h2));
X_h2(:,1) = x0;

% 先计算h/2步长的完整解
for i = 1:length(tspan_h2)-1
    X_h2(:,i+1) = rk4(@(t,x)lorenz(t,x,K,R,r,alpha,mu,gamma), h2, tspan_h2(i), X_h2(:,i));
end

% 计算h步长解，并提取同一时间点的h/2解计算误差
Err = zeros(2, length(tspan));
for i = 1:length(tspan)-1
    X_h(:,i+1) = rk4(@(t,x)lorenz(t,x,K,R,r,alpha,mu,gamma), h, tspan(i), X_h(:,i));
    % 找到h/2时间轴中对应当前时间的索引
    idx = find(tspan_h2 == tspan(i+1));
    if tspan(i+1) <= 200
        Err(1,i+1) = abs(X_h2(1,idx) - X_h(1,i+1));
        Err(2,i+1) = abs(X_h2(2,idx) - X_h(2,i+1));
    end
end

% 筛选t<=200的有效数据，避开t=0
valid_idx = tspan <= 200 & tspan > 0;
Err = Err(:, valid_idx);
t_plot = tspan(valid_idx);

% 绘制误差曲线
plot(log(t_plot), log(Err(1,:)), 'b-', 'LineWidth', 1.5);
hold on;
plot(log(t_plot), log(Err(2,:)), 'r--', 'LineWidth', 1.5);
xlabel('log(时间)');
ylabel('log(误差)');
legend('P的误差', 'Z的误差');
grid on;

内容的提问来源于stack exchange，提问作者frogsforhelp