使用MATLAB ODE15s求解微分方程组遇返回列向量错误及多参数需求

问题背景 #

使用MATLAB的ODE15s求解器求解四阶微分方程组，需要针对7组不同的W和PION参数，存储每组参数下N1、N2、N3、N4的数值，并绘制每组参数下N1、N2、N3、N4随时间的变化曲线。

原代码 #

% Code to solve system of four simultaneous differential equation
clc

N0 = [9.66e17;0;0;0];   % Intital conditions of N1, N2, N3, N4 (9.66e17,0,0,0)
[T, N] = ode15s(@spl_eq,[0 126e-15], N0);  % Intital conditions of time from 0 to 126e-15
plot(T,N(:,4),T,N(:,3),T,N(:,2),T,N(:,1))   % Plotting of N1 N2 N3 N4 vs Time

function dN = spl_eq( t,N )

t = 126e-15;

dN = zeros(4,7);    % Blank array to store N1 N2 N3 N4 from solution for each values of seven values of W and PION
Q = 10e10;      % Constant value
A = 10e8;       % Constant value
Ppd = 10e10;    % Constant value

% W and PION are not constant and they have seven values as row vectors given below:

W = [0.2716 0.4243 0.6110 0.8317 1.0862 1.3748 1.6972]*1e9;
PION = [0.7072 0.8874 1.0254 1.2375 1.4143 1.5911 1.7679]*1e11;


for j = 1:length(W)
    dN(1) = (-W(j)*N(1)) + ((Q+A)*N(2)) + ((W(j)+Q+A)*N(3)) + 0;            % first equation 
    dN(2) = 0 + (-(Q+A)*N(2)) + ((Q+A)*N(3)) + 0;                             % second equation 
    dN(3) = (W(j)*N(1)) + 0 - (W(j)+(2*Q)+(2*A)+ Ppd + PION(j))*N(3) + 0;     % third equation 
    dN(4) = 0 + 0 + PION(j)*N(3) + 0;                                           % fourth equation 
end

end

错误信息 #

Error using odearguments
SPL_EQ must return a column vector.
Error in ode15s (line 148)
    odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in test_soln_rate (line 6)
[T, N] = ode15s(@spl_eq,[0 126e-15], N0);  % Intital conditions of time from 0 to 126e-15

问题分析与修正方案 #

错误原因 #

1. ODE求解器要求导数函数必须返回列向量，但原代码中dN = zeros(4,7)创建了4行7列的矩阵，不符合格式要求。
2. 试图在单个导数函数内同时处理7组参数，这违背了ODE求解器的工作逻辑——求解器每次仅处理一组状态变量，无法同时求解多组参数的方程组。
3. 导数函数内手动赋值t = 126e-15覆盖了求解器传入的时间参数，会导致求解过程完全失效，因为ODE求解器需要根据不同时间点计算导数。

修正后的代码 #

clc

% 定义7组参数
W = [0.2716 0.4243 0.6110 0.8317 1.0862 1.3748 1.6972]*1e9;
PION = [0.7072 0.8874 1.0254 1.2375 1.4143 1.5911 1.7679]*1e11;
N0 = [9.66e17;0;0;0];   % 初始条件
tspan = [0 126e-15];    % 时间范围
Q = 10e10;              % 常量
A = 10e8;               % 常量
Ppd = 10e10;            % 常量

% 预分配存储结果的数组
all_T = cell(1,7);
all_N = cell(1,7);

% 循环求解每组参数的方程组
for j = 1:length(W)
    % 使用匿名函数传递当前组的W和PION参数给导数函数
    [T, N] = ode15s(@(t,N) spl_eq(t,N,W(j),PION(j),Q,A,Ppd), tspan, N0);
    all_T{j} = T;
    all_N{j} = N;
    
    % 绘制当前组的曲线
    figure(j)
    plot(T, N(:,1), 'r-', T, N(:,2), 'g-', T, N(:,3), 'b-', T, N(:,4), 'm-')
    legend('N1','N2','N3','N4')
    xlabel('时间')
    ylabel('数值')
    title(['参数组',num2str(j),': W=',num2str(W(j)),', PION=',num2str(PION(j))])
end

% 导数函数，接收当前时间、状态变量和一组参数
function dN = spl_eq(t,N,W_val,PION_val,Q,A,Ppd)
    dN = zeros(4,1);  % 创建4行1列的列向量，符合ODE求解器要求
    dN(1) = (-W_val*N(1)) + ((Q+A)*N(2)) + ((W_val+Q+A)*N(3));
    dN(2) = -(Q+A)*N(2) + ((Q+A)*N(3));
    dN(3) = W_val*N(1) - (W_val + 2*Q + 2*A + Ppd + PION_val)*N(3);
    dN(4) = PION_val*N(3);
end

修正说明 #

• 将多组参数的求解逻辑移至主程序循环中，每次调用ODE求解器处理一组参数。
• 使用匿名函数@(t,N) spl_eq(t,N,W(j),PION(j),Q,A,Ppd)将当前组的参数传递给导数函数。
• 导数函数内创建zeros(4,1)的列向量返回，严格符合ODE求解器的格式要求。
• 移除导数函数中手动赋值t的错误代码，保留求解器传入的时间参数。
• 预分配cell数组存储每组参数的求解结果，方便后续分析调用。
• 每组参数求解完成后单独绘制曲线，并添加图例和标题区分不同参数组。

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