[Scilab-users] ODE : Second order coupled differential equations

[Stéphane Mottelet]

Hello,

Try to transpose u0.

S.

Le 06/04/2020 à 11:41, Cruz a écrit :
> Hello everybody,
>
> I have a system of 2 differential equations (Eq.1 and Eq.2) coupled of order
> 2 as shown in the attached image. Normally, I can easily find the solution
> by transforming this system onto a system of order 1 by changing variables.
>
> <https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/mailinglists.scilab.org/file/t498030/ode.jpg>
>
> The variables theta and r are functions of time. I use the ODE solver
> integrated in Scilab.
>
> Here is my code :
>
> ------------------------------------------------------------------
>
> t0 =  0; // initial time (s)
> tmax = 30; // solving max time (s)
> tinc = 0.005; // time increment (s)
>
> t = t0:tinc:tmax;
>
> // Differential equations system
>
> function du = equadiff(t,u)
>
>    // d for derivative with respect to time t
>    // u(1) = theta
>    // u(2) = d(theta)/dt
>    // u(3) = r
>    // u(4) = d(r)/dt
>    // du(1) = d(theta)/dt
>    // du(2) = d²(theta)dt²
>    // du(3) = d(r)/dt
>    // du(4) = d²(r)/dt²
>
>    F1=1300;
>    F2=1000;
>    l=13;
>    k=68000;
>    menv = 112.14;
>    m0 = 56.07;
>    g=9.81;
>
>    m = (menv - m0) * (1-abs(cos(u(1)))) + m0;
>
>    du(1) = u(2);
>    du(2) = (F2*sin(u(1)) + (F1-m*g).*cos(u(1)) - 2*m.*u(4).*u(2)) ./
> (m.*u(3));
>    du(3) = u(4);
>    du(4) = ((F1-m*g).*sin(u(1)) - F2*cos(u(1)) - k*(u(3)-l) +
> m.*u(3).*u(2).*u(2)) ./ m;
>
> endfunction
>
> // Integration
>
> u0 = [0,0,13,0];
> [u] = ode(u0,t0,t,equadiff);
>
> ------------------------------------------------------------------
>
> The problem is that the output of the solver, [u], should be a matrix of 4
> lines and 6001 columns:
>
> 4 lines:
>
> u (1) = theta
> u (2) = d (theta) / dt
> u (3) = r
> u (4) = d (r) / dt
>
> 6001 columns:
>
> tmax / tinc = 6000 time step + the initial instant t0 = 6001
>
> At the end of the calculation I get
>
> size (u)
>
> ans =
>
> 1. 24004.
>
> instead of
>
> 4. 6001.
>
> It seems to me to proceed in the standard way for this type of problem.
> Someone have any idea why it doesn't work?
>
> Thanks in advance to all !
>
>
>
>
>
>
>
>
>
>
> --
> Sent from: https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
> _______________________________________________
> users mailing list
> users at lists.scilab.org
> https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users

-- 
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
http://www.utc.fr/~mottelet