如何求解耦合非线性一阶微分方程并在MATLAB中绘图？

阿华AIGC实验室

2026-5-21

求解耦合微分方程并绘制x-t、y-t曲线（MATLAB实现）

看起来你遇到的是一组特殊的耦合微分方程——方程两边互相包含对方的导数，直接用常规ODE求解器会有问题。我帮你梳理下处理思路，再给出具体的MATLAB实现方案。

步骤1：整理方程形式

你的两个方程：

dx/dt = (-a1sin(y) + a2 + a3sin(y-x)) / ((dy/dt)a4cos(y-x))
dy/dt = (a1sin(x) -a5 + a6x + a7*sin(y-x)) / ((dx/dt)a8cos(y-x))

首先把两个方程交叉相乘，消去分母中的导数项，会得到一个代数约束关系：

(-a1*sin(y) + a2 + a3*sin(y-x)) / a4 = (a1*sin(x) -a5 + a6*x + a7*sin(y-x)) / a8

这说明x(t)和y(t)并非完全独立，它们始终满足这个代数等式。结合初始条件（比如x(t₀)=x₀，y(t₀)=y₀，且初始值必须满足上述约束），我们可以把问题转化为微分代数方程组（DAE），或者进一步简化为单变量ODE。

不过这里先提一种可能性：如果方程中的(dy/dt)和(dx/dt)是笔误（比如原本分母只有a4*cos(y-x)和a8*cos(y-x)），那方程会变成标准的一阶耦合ODE组，求解会简单很多。我先把这种情况的方案列出来，再处理原方程的情况。

方案1：假设方程存在笔误（标准ODE组）

如果原方程实际是：

dx/dt = (-a1*sin(y) + a2 + a3*sin(y-x)) / (a4*cos(y-x))
dy/dt = (a1*sin(x) -a5 + a6*x + a7*sin(y-x)) / (a8*cos(y-x))

这是典型的一阶ODE组，可以直接用MATLAB的ode45求解，步骤如下：

代码实现

定义ODE函数：

function dydt = myODE(t, z, a1, a2, a3, a4, a5, a6, a7, a8)
    x = z(1);
    y = z(2);
    % 计算dx/dt
    dxdt = (-a1*sin(y) + a2 + a3*sin(y-x)) / (a4*cos(y-x));
    % 计算dy/dt
    dydt = (a1*sin(x) - a5 + a6*x + a7*sin(y-x)) / (a8*cos(y-x));
    dydt = [dxdt; dydt];
end

设置参数、初始条件和求解时间：

% 设置参数a1-a8（替换为你的实际值）
a1 = 1; a2 = 0.5; a3 = 0.3; a4 = 2;
a5 = 0.2; a6 = 0.1; a7 = 0.4; a8 = 1.5;

% 初始条件t0=0时的x0和y0
t0 = 0;
x0 = pi/4;
y0 = pi/6;
z0 = [x0; y0];

% 求解时间范围，比如从0到10
tspan = [0 10];

% 调用ode45求解
[t, z] = ode45(@(t,z) myODE(t,z,a1,a2,a3,a4,a5,a6,a7,a8), tspan, z0);

% 提取x(t)和y(t)
x = z(:,1);
y = z(:,2);

绘制x-t和y-t曲线：

figure;
subplot(2,1,1);
plot(t, x, 'b-', 'LineWidth', 1.5);
xlabel('t'); ylabel('x(t)');
title('x-t曲线');
grid on;

subplot(2,1,2);
plot(t, y, 'r-', 'LineWidth', 1.5);
xlabel('t'); ylabel('y(t)');
title('y-t曲线');
grid on;

方案2：按照原方程处理（微分代数方程组DAE）

如果原方程的形式完全正确，我们需要处理DAE问题（包含微分方程和代数约束）。MATLAB的ode15i专门用于求解隐式ODE（包括DAE），步骤如下：

步骤说明

首先把原方程转化为隐式形式F(t, z, zdot) = 0，其中z = [x; y]，zdot = [dxdt; dydt]：

zdot(1)*zdot(2)*a4*cos(z(2)-z(1)) + a1*sin(z(2)) - a2 - a3*sin(z(2)-z(1)) = 0
zdot(1)*zdot(2)*a8*cos(z(2)-z(1)) - a1*sin(z(1)) + a5 - a6*z(1) - a7*sin(z(2)-z(1)) = 0

初始条件需要满足：

z0 = [x0; y0]（必须满足代数约束）
需要提供初始导数的猜测值zdot0（可以通过代数约束和其中一个方程估算）

代码实现

定义隐式ODE函数：

function F = myDAE(t, z, zdot, a1, a2, a3, a4, a5, a6, a7, a8)
    x = z(1);
    y = z(2);
    dxdt = zdot(1);
    dydt = zdot(2);
    
    % 第一个方程的隐式形式
    F1 = dxdt*dydt*a4*cos(y-x) + a1*sin(y) - a2 - a3*sin(y-x);
    % 第二个方程的隐式形式
    F2 = dxdt*dydt*a8*cos(y-x) - a1*sin(x) + a5 - a6*x - a7*sin(y-x);
    
    F = [F1; F2];
end

设置参数、初始条件并求解：

% 设置参数（替换为你的实际值）
a1 = 1; a2 = 0.5; a3 = 0.3; a4 = 2;
a5 = 0.2; a6 = 0.1; a7 = 0.4; a8 = 1.5;

% 初始条件t0=0时的x0和y0，必须满足代数约束
t0 = 0;
x0 = pi/4;
% 用fsolve求解满足约束的y0（初始猜测值设为pi/6）
fun_y = @(y) (-a1*sin(y)+a2+a3*sin(y-x0))/a4 - (a1*sin(x0)-a5+a6*x0+a7*sin(y-x0))/a8;
y0 = fsolve(fun_y, pi/6);

z0 = [x0; y0];

% 求解初始导数的猜测值（初始猜测值设为[0.1; 0.1]）
fun_zdot = @(zdot) myDAE(t0, z0, zdot, a1,a2,a3,a4,a5,a6,a7,a8);
zdot0 = fsolve(fun_zdot, [0.1; 0.1]);

% 求解时间范围
tspan = [0 10];

% 调用ode15i求解
opts = odeset('RelTol',1e-6,'AbsTol',1e-8);
[t, z] = ode15i(@(t,z,zdot) myDAE(t,z,zdot,a1,a2,a3,a4,a5,a6,a7,a8), tspan, z0, zdot0, opts);

% 提取x(t)和y(t)
x = z(:,1);
y = z(:,2);

绘制曲线（和方案1的绘图代码完全相同）：

figure;
subplot(2,1,1);
plot(t, x, 'b-', 'LineWidth', 1.5);
xlabel('t'); ylabel('x(t)');
title('x-t曲线');
grid on;

subplot(2,1,2);
plot(t, y, 'r-', 'LineWidth', 1.5);
xlabel('t'); ylabel('y(t)');
title('y-t曲线');
grid on;