寻求MATLAB求解两个联立常微分方程的技术指导

问题明确 #

给定参数：

DKS = 1.96e-9;
DLS = 1.03e-9;
JT = 0.000506707;
JK = 0.000268163;
% 自变量x范围：0到0.0005 m
% 初始条件：x=0时，K(0)=125，L(0)=125

未知量为K(x)和L(x)，修正原方程符号笔误后，对应的微分方程为：

ODE1：$-DKS \cdot \frac{dK}{dx} \cdot \frac{2K+L}{K+L} - DKS \cdot \frac{K}{K+L} \cdot \frac{dL}{dx} - JK = 0$
ODE2：$-DLS \cdot \frac{dL}{dx} \cdot \frac{K+2L}{K+L} + DLS \cdot \frac{L}{K+L} \cdot \frac{d^2K}{dx2} + (JT-JK) = 0$

关键转换：二阶ODE降为一阶方程组 #

ode45仅支持求解一阶常微分方程组，因此需要将含二阶导数的ODE2转换为一阶形式：

设变量替换：

• $y_1 = K(x)$
• $y_2 = \frac{dK}{dx}$
• $y_3 = L(x)$
• $y_4 = \frac{dL}{dx}$

推导一阶导数表达式 #

1. 从ODE1解出$\frac{dL}{dx}$：
   $$\frac{dL}{dx} = -\frac{JK \cdot (K+L) + DKS \cdot \frac{dK}{dx} \cdot (2K+L)}{DKS \cdot K}$$
   对应变量替换：$y_4 = -\frac{JK \cdot (y_1+y_3) + DKS \cdot y_2 \cdot (2y_1+y_3)}{DKS \cdot y_1}$

2. 从ODE2解出$\frac{d^2K}{dx2}$（即$\frac{dy_2}{dx}$）：
   $$\frac{d^2K}{dx2} = \frac{DLS \cdot \frac{dL}{dx} \cdot (K+2L) - (JT-JK) \cdot (K+L)}{DLS \cdot L}$$
   代入$y_4 = \frac{dL}{dx}$后：
   $$\frac{dy_2}{dx} = \frac{DLS \cdot y_4 \cdot (y_1+2y_3) - (JT-JK) \cdot (y_1+y_3)}{DLS \cdot y_3}$$

3. 剩余一阶导数关系：
   $$\frac{dy_1}{dx} = y_2, \quad \frac{dy_3}{dx} = y_4$$

MATLAB代码实现 #

1. 定义微分方程组函数 #

创建名为ode_system.m的函数文件：

function dydx = ode_system(x, y, DKS, DLS, JT, JK)
    % 变量映射
    y1 = y(1); % K(x)
    y2 = y(2); % dK/dx
    y3 = y(3); % L(x)
    
    % 计算y4（dL/dx），从ODE1推导
    numerator_y4 = JK * (y1 + y3) + DKS * y2 * (2*y1 + y3);
    y4 = -numerator_y4 / (DKS * y1);
    
    % 计算dy2/dx（d²K/dx²），从ODE2推导
    numerator_dy2 = DLS * y4 * (y1 + 2*y3) - (JT - JK) * (y1 + y3);
    dy2dx = numerator_dy2 / (DLS * y3);
    
    % 组装一阶方程组的导数向量
    dydx = [
        y2;          % dy1/dx = dK/dx = y2
        dy2dx;       % dy2/dx = d²K/dx²
        y4;          % dy3/dx = dL/dx = y4
        0;           % y4的导数由y1/y2/y3的导数自动推导，无需显式定义
    ];
end

2. 主脚本求解与绘图 #

创建主脚本文件（例如solve_ode.m）：

% 定义参数
DKS = 1.96e-9;
DLS = 1.03e-9;
JT = 0.000506707;
JK = 0.000268163;

% 自变量x的求解范围
xspan = [0, 0.0005];

% 初始条件：y(1)=K(0)=125, y(2)=dK/dx(0), y(3)=L(0)=125, y(4)=dL/dx(0)
% 原问题未给出一阶导数初始值，这里假设为零通量边界（导数为0），可根据实际情况修改
y0 = [125; 0; 125; 0];

% 使用ode45求解，传入额外参数
[x, y] = ode45(@(x,y) ode_system(x,y,DKS,DLS,JT,JK), xspan, y0);

% 绘图：K(x)和L(x)随x的变化
figure;
plot(x, y(:,1), 'b-', 'LineWidth', 1.5);
hold on;
plot(x, y(:,3), 'r--', 'LineWidth', 1.5);
xlabel('x (m)');
ylabel('Concentration (moles/m^3)');
legend('K(x)', 'L(x)');
grid on;
title('K and L Concentration Distribution');

注意事项 #

• 初始条件补充：原问题只给出了K(0)和L(0)，但转换后的一阶方程组需要dK/dx(0)和dL/dx(0)的初始值。代码中假设为零通量（导数为0），如果实际问题有其他边界条件，请修改y0中的对应值，或通过x=0处的ODE1推导初始导数。
• 方程符号检查：原代码中的CLB、CKB疑似应为未知量L(x)和K(x)，若实际符号不同，请调整代码中的变量映射。
• 数值稳定性：由于参数数量级差异较大（如DKS是1e-9，JT是1e-4），求解时若出现警告，可尝试使用ode15s等刚性求解器替代ode45。

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