在MATLAB中使用ode45求解含4D矩阵的二阶常微分方程

问题背景 #

给定二阶矩阵微分方程：
B1*phi + B2*(dphi/ds) - B3*(d²phi/ds²) = 0
其中B1、B2、B3为4×4矩阵，初始条件为s=0时，phi=0（4维零向量），dphi/ds=0（4维零向量）。

用户提供的代码 #

% Define the anonymous function for the equation
equation = @(s, phi_dphi) [phi_dphi(2); -B1.*phi_dphi(1) + B2.*phi_dphi(2) - B3.*(phi_dphi(2).^2)];

% Define the initial conditions
phi0 = [0;0;0;0]; % Initial condition for phi
dphi_ds0 = [0;0;0;0]; % Initial condition for d(phi)/ds

% Set the options for ode45 (optional)
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);

% Solve the equation using ode45
[t, phi_dphi] = ode45(equation, [S_0, S_L], [phi0; dphi_ds0], options);

% Extract the solution for phi and d(phi)/ds
phi = phi_dphi(:, 1);
dphi_ds = phi_dphi(:, 2);

% Display the results
disp('Solution:')
disp('---------')
disp('s phi d(phi)/ds')
disp([t, phi, dphi_ds]);

代码存在的问题 #

1. 一阶方程组转化错误
   原方程整理为d²phi/ds² = B3⁻¹*(B1*phi + B2*dphi/ds)（需保证B3可逆），但用户代码里的右侧表达式完全不符合推导逻辑，错误使用了元素乘.*和平方项，完全偏离原方程。

2. 状态变量维度处理错误
   phi是4维向量，因此状态变量phi_dphi应该是8维（前4维为phi，后4维为dphi/ds），但用户代码里仅取了单个元素phi_dphi(2)，无法表示4维向量的导数。

3. 结果提取错误
   求解后phi_dphi的每一行是对应s处的8维状态，phi应取前4列，dphi/ds取后4列，用户仅提取第1、2列，丢失了大部分维度信息。

4. 未定义变量
   代码中S_0和S_L未提前赋值，运行会直接报错。

修正后的代码 #

% 提前定义4×4矩阵B1、B2、B3（示例，需替换为实际矩阵）
B1 = rand(4);
B2 = rand(4);
B3 = rand(4);
% 确保B3可逆，若不可逆需特殊处理
if det(B3) < 1e-10
    error('B3矩阵不可逆，无法直接求逆');
end
B3_inv = inv(B3);

% 定义一阶方程组：y = [phi; dphi/ds]，则dy/ds = [dphi/ds; d²phi/ds²]
equation = @(s, y) [y(5:8); B3_inv*(B1*y(1:4) + B2*y(5:8))];

% 初始条件
phi0 = zeros(4,1);
dphi_ds0 = zeros(4,1);
y0 = [phi0; dphi_ds0];

% 求解区间（需替换为实际的起始和终止s值）
S_0 = 0;
S_L = 10;

% ODE求解选项
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);

% 求解ODE
[s_vals, y_sol] = ode45(equation, [S_0, S_L], y0, options);

% 提取结果
phi_sol = y_sol(:, 1:4); % 每一行对应一个s的4维phi
dphi_ds_sol = y_sol(:, 5:8); % 每一行对应一个s的4维dphi/ds

% 显示部分结果（因维度较高，仅展示前5行）
disp('s值 | phi前2维 | dphi/ds前2维');
disp('-----------------------------');
disp([s_vals(1:5), phi_sol(1:5,1:2), dphi_ds_sol(1:5,1:2)]);

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