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);
3 个评论
Jan
2019-5-9
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
2019-5-9
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
2019-5-9
编辑：Davide Figoli
2019-5-9
回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息
2019-5-8
2019-5-9
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
