Solving second order differential system using ode45

Hi there every-one,

I'm tryng to resolve a second order differential system using ode45.

I've tried to write the system in this way and i resolved it using ode45 but i'm not conviced that i've done it correctly.

The function torque is a function that gived a time return the values of torque gived by a motor.

I would like that that ode45 will return every ddth, dth and th solution of the system.

function yd = dinamic_flex(t, y)
    
    global J1 J2 Jmr J3 J4 r1 r2 r3 r4 ...
        kgiu kcin12 kcin34 
    
   Mm= torque(t);

    th2 = y(1); th1 = y(2); thm = y(3); th3 = y(4); th4 = y(5);   
    dth2 = y(6); dth1 = y(7); dthm = y(8); dth3 = y(9); dth4 = y(10);

    ddth2 = (1 / J2) * (r2 * kcin12 * ((r1 * th1) - (r2 * th2)) + 0 * (dth1 - dth2));
    ddth1 = (1 / J1) * ((thm - th1) * 2*kgiu + kcin12 * r1 *((r2 * th2) - (r1 * th1)) + 0 * (dthm - dth1) + 0 * (dth2 - dth1));
    ddthm = (1 / Jmr) * ((th1 - thm) * 2*kgiu + (th3 - thm) * 2*kgiu + 0 * (dth1 - dthm) + 0 * (dth3 - dthm) + Mm);
    ddth3 = (1 / J3) * ((thm - th3) * 2*kgiu + r3 * kcin34 * ((r4 * th4) - (r3 * th3)) + 0 * (dthm - dth3) + 0 * (dth4 - dth3));
    ddth4 = (1 / J4) * (r4 * kcin34 * ((th3 * r3) - (th4 * r4)) - 0 + 0 * (dth3 - dth4));
    
    yd(1) = dth2; yd(2) = dth1; yd(3) = dthm; yd(4) = dth3; yd(5) = dth4;   
    yd(6) = ddth2; yd(7) = ddth1; yd(8) = ddthm; yd(9) = ddth3; yd(10) = ddth4;
    
    yd = yd';
    
end

 
yzero = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0];
op = 1e-9;
options = odeset('RelTol', op);

[t, y_f] = ode45(@dinamica_flex, t_analysis, yzero, options);

[Jan]

All we see is the current code. This code computes, what this code computes. We do not have any chance to estimate, if the code is correct or not, because there is no explanation, what it should compute. So what exactly is your question?

[David Wilson]

I'm making a few guesses here. I re-wrote your internal ODE as the following:

function yd = dinamic_flex(t, y, ... 
              J1, J2, Jmr, J3, J4, r1, r2, r3, r4, ...
              kgiu, kcin12, kcin34)
 
   %Mm= torque(t); % don know exaclty what this is so ... 
   Mm = sin(t); % fake it 
   
    th2 = y(1); th1 = y(2); thm = y(3); th3 = y(4); th4 = y(5);   
    dth2 = y(6); dth1 = y(7); dthm = y(8); dth3 = y(9); dth4 = y(10);
    ddth2 = (1 / J2) * (r2 * kcin12 * ((r1 * th1) - (r2 * th2)) + 0 * (dth1 - dth2));
    ddth1 = (1 / J1) * ((thm - th1) * 2*kgiu + kcin12 * r1 *((r2 * th2) - (r1 * th1)) + 0 * (dthm - dth1) + 0 * (dth2 - dth1));
    ddthm = (1 / Jmr) * ((th1 - thm) * 2*kgiu + (th3 - thm) * 2*kgiu + 0 * (dth1 - dthm) + 0 * (dth3 - dthm) + Mm);
    ddth3 = (1 / J3) * ((thm - th3) * 2*kgiu + r3 * kcin34 * ((r4 * th4) - (r3 * th3)) + 0 * (dthm - dth3) + 0 * (dth4 - dth3));
    ddth4 = (1 / J4) * (r4 * kcin34 * ((th3 * r3) - (th4 * r4)) - 0 + 0 * (dth3 - dth4));
    
    yd(1) = dth2; yd(2) = dth1; yd(3) = dthm; yd(4) = dth3; yd(5) = dth4;   
    yd(6) = ddth2; yd(7) = ddth1; yd(8) = ddthm; yd(9) = ddth3; yd(10) = ddth4;
    
    yd = yd(:); % force column 
    
end

I removed the globals (never a good idea anyway), and substituted your torque function for something that I can compute. You will of course have to change that back again.

Now since you didn't tell us what the constants were, I made them up. I then called ode45 to solve the system starting from the origin, and simulated for 10 time units.

% Set some constants ; 
J1 =1; J2=2; Jmr=3; J3=4; J4=5; r1=6; r2=7; r3=8; r4=9;
kgiu=10; kcin12=11; kcin34=13; 

yzero = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]';
tol = 1e-9;
options = odeset('RelTol', tol);
t_analysis = [0,10]; % guessing here ... 

[t, y_f] = ode45(@(t,y) dinamic_flex(t, y, ... 
              J1, J2, Jmr, J3, J4, r1, r2, r3, r4, kgiu, kcin12, kcin34), ... 
              t_analysis, yzero, options); % spell function correctly here


plot(t, y_f) 

I get some trends, and I guess it is up to you to see if they are sensible. By the eay, you spelt the name of the function, dinamic(a)_flex, two different ways.

[Davide Figoli]

Thank you very much for the advice regarding the globals.

I apologize maybe i didn't have explained my-self correctelly.

th2, th1,thm,th3,th4 represent angles, so dth2, dth1,dthm,dth3,dth4 are angular velocity and ddth2, ddth1,ddthm,ddth3,ddth4 are angular accelleration of a mechanicle system. I use:

Mm= torque(t)

To give to the cranck of the engine a torque dipending by the time. In fact Mm is used in the equation regarding the engine ( ddthm,dthm, thm the letter m stands for motor)

 ddthm = (1 / Jmr) * ((th1 - thm) * 2*kgiu + (th3 - thm) * 2*kgiu + 0 * (dth1 - dthm) + 0 * (dth3 - dthm) + Mm);

My question is is possible the restitution of all the values ( th2, th1,thm,th3,th4,dth2, dth1,dthm,dth3,dth4,ddth2, ddth1,ddthm,ddth3,ddth4) that solves the system by the ode45?

[t, y_f] = ode45(@dinamica_flex, t_analysis, yzero, options);

In other words i would like to have y_f as a (length(t),15) vector where the first 5 columns represent my th2, th1,thm,th3,th4 the colums form 6 to 10 represent dth2, dth1,dthm,dth3,dth4 and the colums from 10 to 15 represent ddth2, ddth1,ddthm,ddth3,ddth4.

Thank you very much.