Solving second order differential system using ode45
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
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);
3 Kommentare
David Wilson
am 9 Mai 2019
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.
Antworten (0)
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!