How to solve a set of ODEs in Scilab?

Question

I have following set of ODEs which I would like to numerically solve in Scilab

I have successfully written the function for evaluating the right hand side of the first equation in other words I am able to solve the first equation

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

My problem is that I don't know how to write similar function for evaluation of the right hand side of the second equation because it depends on solution of the first equation. Can anybody give me an advice how to do that?

Possible solution:

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));

What is x_tilde with respect to x? Could you name it `y` instead? A and B look constant, aren't they? If so, you can define them outside SystemModel(), hence only once for all. Are An, L, Cn, ΔA, ΔB and ΔC constant wrt t? And are they scalar, or matrices? — S. Gougeon, Feb 12 '23 at 22:11
@S.Gougeon thank you for your reaction. As far as your questions. x_tilde is independent in respect to x. All the matrices (A, B, An, L, Cn, delta_A, delta_C, delta_B) contains constant elements. — Steve, Feb 12 '23 at 22:15
To me, you should first completely expand your system of ODEs: replace x (4x1) with [x1;x2;x3;x3], dx/dt with [dx1/dt; dx2/dt; dx3/dt; dx4/dt], etc, and build the explicit system with only scalar unknowns. — S. Gougeon, Feb 12 '23 at 22:35
... and simplified the notations. for instance: D = An - L.Cn ; E = ΔA - L.ΔC ; and express things with D and E predefined coordinates, not in a matricial way. — S. Gougeon, Feb 12 '23 at 22:44

S. Gougeon · Accepted Answer · 2023-02-12T23:30:28.647

1

Let y = x_tilde. Let assume that it is a 3x1 vector (we can't guess its size with your current presentation).

Build the column X = [x1 x2 x3 x4 y1 y2 y3].' (big X)
Express the column (dX/dt) according to X coordinates and t
Convert the system built in 2) into a Scilab function X_dot = Xder(t, X)
Build the initial state vector Xinit = [x1(t_init); x2(t_init); .. y3(t_init)]
Define the vector t of times to which you want the values of X. They all likely have to be ≥ t_init, and must be strictly increasing.
call X = ode(Xinit, t_init, t, Xder)

X(:,i) should then be the values of X components for each t(i) date.

You can "back-split" big X into x = X(1:4,:) and x_tilde = X(5:$,:).

edited Feb 12 '23 at 23:30

answered Feb 12 '23 at 23:14

S. Gougeon

791
3
16

Thank you very much for your answer. The approach you have proposed above lead me to the idea I have attempted to describe as a part of my question. The motivation for the approach I have describe above is that for me it's more convenient to work with the ODEs in blocks. What do you think about this "improvement"? – Steve Feb 13 '23 at 08:47
The main thing to use `ode()` is to build a unique vectorized unknow by concatenating x and x_tilde. Within `RightHandSide()`, indeed you may split `X` and then formulate the system as you do, in a matricial way. – S. Gougeon Feb 13 '23 at 10:07
1

But, still: all your heavy calculations of constants from `A=..` to `dC=..` should be outside RightHandSide(), since they do not depend neither on X nor on t,. `RightHandSide()` is called by ode() a lot of times. Computing all these coefficients 1e6 times while they are constant is very time-consuming and meaningless. You can leave them in the calling environment, or include them in the calling `list(RightHandSide, A, ..)`. By the way, putting u as a RightHandSide's argument is useless (but not really harmful), because u() is a public function, as soon as you define it before calling ode() – S. Gougeon Feb 13 '23 at 10:08

How to solve a set of ODEs in Scilab?

1 Answers1