MATLAB中ode45指定Mass矩阵时odeset参数设置报错问询

问题场景 #

我有一个对应非线性系统的ODE方程，编写MATLAB求解代码后，未添加options参数时ode45可正常运行，但调用ode45时传入包含Mass矩阵与容差设置的odeset参数，程序出现报错。需要找到正确设置odeset参数适配ode45的方法。

原报错代码示例 #

主求解代码

% 时间区间
tspan = [0 10];
% 初始条件
y0 = [1; 0];
% 构建Mass矩阵
M = [1 0; 0 2];
% 设置options
options = odeset('Mass', M, 'RelTol', 1e-6, 'AbsTol', 1e-8);
% 调用ode45
[t, y] = ode45(@myOdeFun, tspan, y0, options);

原myOdeFun函数

function dydt = myOdeFun(t, y)
    % 非线性ODE定义
    dydt = [y(2); -sin(y(1))];
end

正确设置方案 #

核心原理 #

当通过odeset设置Mass参数时，ode45会将方程视为**M(t,y)·y’ = f(t,y)**的形式（默认情况下等价于M为单位矩阵，此时y’ = f(t,y)）。因此必须调整odefun的输出，使其对应f(t,y)而非直接返回y’。

方案1：直接调整odefun输出适配固定Mass矩阵 #

修改myOdeFun，让其返回满足M·y’ = f(t,y)的f(t,y)：

function f = myOdeFun(t, y)
    % 原方程y' = [y(2); -sin(y(1))]，结合M=[1 0;0 2]，推导f = M*y'
    f = [y(2); -2*sin(y(1))];
end

主程序无需修改，直接运行即可正常调用ode45。

方案2：将Mass矩阵定义为动态函数传入 #

如果Mass矩阵随时间t或状态y变化，可将其定义为单独函数，通过odeset传入：

% 主程序修改options设置
options = odeset('Mass', @massMatrixFun, 'RelTol', 1e-6, 'AbsTol', 1e-8);

% 定义动态Mass矩阵函数
function M = massMatrixFun(t, y)
    % 这里可根据t或y动态生成矩阵，示例为固定矩阵
    M = [1 0; 0 2];
end

同样需要确保myOdeFun返回的是f(t,y)，满足M(t,y)·y’ = f(t,y)。

验证代码 #

修改后运行以下代码可验证结果：

tspan = [0 10];
y0 = [1; 0];
M = [1 0; 0 2];
options = odeset('Mass', M, 'RelTol', 1e-6, 'AbsTol', 1e-8);
[t, y] = ode45(@myOdeFun, tspan, y0, options);

% 绘图查看结果
plot(t, y(:,1), 'r-', t, y(:,2), 'b--');
legend('y_1', 'y_2');
xlabel('时间');
ylabel('状态量');

报错原因 #

原代码中odefun直接返回y’，对应默认的y’ = f(t,y)形式，但传入Mass参数后，ode45期望方程遵循M·y’ = f(t,y)的格式，两者输出不匹配导致程序报错。

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