如何在Scilab中求解存在依赖关系的常微分方程组？

需要在Scilab中数值求解一组耦合常微分方程组，已实现第一个方程的右侧计算函数并能独立求解，但第二个方程的求解依赖第一个方程的解，不清楚如何编写对应的右侧函数。以下是现有代码及可行解决方案：

现有代码 #

function ut = u(t)
    ut = [Vm*cos(2*%pi*fs*t); Vm*sin(2*%pi*fs*t)];
endfunction

function dxdt = SystemModel(t, x, u)
    A = [(-RS*(LL + LM)^2 - RR*LM^2)/(LL*LM*(LL + LM)),                                                       0.0, RR/(LL*(LL + LM)),             pp*wm/LL;
                                                   0.0,             (-RS*(LL + LM)^2 - RR*LM^2)/(LL*LM*(LL + LM)),         -pp*wm/LL,    RR/(LL*(LL + LM));
                                     (LM*RR)/(LL + LM),                                                       0.0,     -RR/(LL + LM),               -pp*wm;
                                                   0.0,                                         (LM*RR)/(LL + LM),             pp*wm,        -RR/(LL + LM)];
                                                         
    B = [(LL + LM)/(LL*LM), 0.0; 0.0, (LL + LM)/(LL*LM); 0.0, 0.0; 0.0, 0.0];
    
    dxdt = A*x + B*u(t);
endfunction

可行解决方案代码 #

x0 = zeros(4, 1);
xtilde0 = zeros(4, 1);
X0 = [x0; xtilde0];
t0 = 0;
dt = 0.001;
t = 0:dt:1;

function ut = u(t)
    ut = [Vm*cos(2*%pi*fs*t); Vm*sin(2*%pi*fs*t)];
endfunction

function dXdt = RightHandSide(t, X, u)

    x      = X(1:4);
    xtilde = X(5:8);

    // dx/dt = A*x + B*u 
    A = [(-RS*(LL + LM)^2 - RR*LM^2)/(LL*LM*(LL + LM)),                                                       0.0, RR/(LL*(LL + LM)),             pp*wm/LL;
                                                   0.0,             (-RS*(LL + LM)^2 - RR*LM^2)/(LL*LM*(LL + LM)),         -pp*wm/LL,    RR/(LL*(LL + LM));
                                     (LM*RR)/(LL + LM),                                                       0.0,     -RR/(LL + LM),               -pp*wm;
                                                   0.0,                                         (LM*RR)/(LL + LM),             pp*wm,        -RR/(LL + LM)];
                                                         
    B = [(LL + LM)/(LL*LM), 0.0; 0.0, (LL + LM)/(LL*LM); 0.0, 0.0; 0.0, 0.0];
    
    // dxtilde/dt = (An - L*Cn)*xtilde + (dA - L*dC)*x + dB*u
    An = [(-RSn*(LLn + LMn)^2 - RRn*LMn^2)/(LLn*LMn*(LLn + LMn)),                                                     0.0, RRn/(LLn*(LLn + LMn)),             pp*wm/LLn;
                                                             0.0,  (-RSn*(LLn + LMn)^2 - RRn*LMn^2)/(LLn*LMn*(LLn + LMn)),            -pp*wm/LLn, RRn/(LLn*(LLn + LMn));
                                           (LMn*RRn)/(LLn + LMn),                                                     0.0,      -RRn/(LLn + LMn),                -pp*wm;
                                                             0.0,                                   (LMn*RRn)/(LLn + LMn),                 pp*wm,     -RRn/(LLn + LMn)];
                                                              
    K = 1.5;
    l1 = (K - 1.0)*((RSn*(LLn + LMn)^2 + RRn*LMn^2)/(LLn*LMn*(LLn + LMn)) + RRn/(LLn + LMn));
    l2 = (K - 1.0)*pp*wm;
    l3 = (K^2 - 1.0)*((RSn*(LLn + LMn)^2 + RRn*LMn^2)/(LMn*(LLn + LMn)) - (LMn*RRn)/(LLn + LMn)) - (K - 1)*((RSn*(LLn + LMn)^2 + RRn*LMn^2)/(LMn*(LLn + LMn)) + (LLn*RRn)/(LLn + LMn));
    l4 = -(K - 1.0)*LLn*wm*pp;
    L = [l1, l2; 
        -l2, l1; 
         l3, l4;
        -l4, l3];
        
    Bn = [(LLn + LMn)/(LLn*LMn), 0.0; 0.0, (LLn + LMn)/(LLn*LMn); 0.0, 0.0; 0.0, 0.0];
                                                              
    Cn = [1.0, 0.0, 0.0, 0.0; 0.0, 1.0, 0.0, 0.0];
    
    A = [(-RS*(LL + LM)^2 - RR*LM^2)/(LL*LM*(LL + LM)),                                                       0.0, RR/(LL*(LL + LM)),             pp*wm/LL;
                                                   0.0,             (-RS*(LL + LM)^2 - RR*LM^2)/(LL*LM*(LL + LM)),         -pp*wm/LL,    RR/(LL*(LL + LM));
                                     (LM*RR)/(LL + LM),                                                       0.0,     -RR/(LL + LM),               -pp*wm;
                                                   0.0,                                         (LM*RR)/(LL + LM),             pp*wm,        -RR/(LL + LM)];
                                                         
    B = [(LL + LM)/(LL*LM), 0.0; 0.0, (LL + LM)/(LL*LM); 0.0, 0.0; 0.0, 0.0];
    
    C = [1.0, 0.0, 0.0, 0.0; 0.0, 1.0, 0.0, 0.0];
    
    dA = An - A;
    dB = Bn - B;
    dC = Cn - C;
    
    dxdt      = A*x + B*u(t);
    dxtildedt = (An - L*Cn)*xtilde + (dA - L*dC)*x + dB*u(t);
    dXdt = [dxdt; dxtildedt];
endfunction

X = ode(X0, t0, t, list(RightHandSide, u));

关键实现思路 #

• 状态向量合并：将两个方程组的状态变量x和xtilde合并为一个大的状态向量X，让Scilab的ode函数可以统一处理整个耦合系统
• 状态拆分：在右侧计算函数中，从合并后的X拆分出各自的子状态向量，用于分别计算每个方程组的导数
• 耦合计算：第二个方程的导数dxtildedt直接调用第一个方程的当前状态x，实现两个方程组的耦合求解
• 统一返回导数：将两个子导数向量合并为大的导数向量dXdt返回，供ode函数完成数值积分

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