Azzera filtri
Azzera filtri

Solve a system of 3 second order equations

2 visualizzazioni (ultimi 30 giorni)
Sam Butler
Sam Butler il 7 Feb 2023
Commentato: Star Strider il 8 Feb 2023
I want to solve the following 3 equations using MATLAB:
The unknowns to the three equations are theta(1), theta(2) and theta(n). Theta(out) is the motor input speed which is defined by a periodic equation however, for the sake of this, can be given some arbitary value. The rest are constants.
I want to solve these over a period of 0 to 5 seconds in the aim to produce torque deflection graphs.
I have searched extensivley for the best way to do it and after multiple trials, I am still struggling (ode45, odefun etc.)
Any help would be greatly appreciated.

Risposte (3)

Alan Stevens
Alan Stevens il 7 Feb 2023
Write your equations as 6 first order equations like so:
where I've assumed you have pre-defined omegaout and thetaout as functions of time.
You will also need initial conditions for the thetas and omegas.
Then use ode45.

Star Strider
Star Strider il 7 Feb 2023
This should get you started —
syms c_1 c_2 k_1 k_2 k_n k_nl J_1 J_2 J_n theta_1(t) theta_2(t) theta_n(t) theta_out(t) t Y
d1theta_1 = diff(theta_1);
d2theta_1 = diff(theta_1,2);
d1theta_2 = diff(theta_2);
d2theta_2 = diff(theta_2,2);
d1theta_n = diff(theta_n);
d2theta_n = diff(theta_n,2);
d1theta_out = diff(theta_out);
DE1 = J_1*d2theta_1 + c_1*(d1theta1 - d1theta_out) + k_1*(theta_1 - theta_out) +
DE2 =
DE3 =
[VF, Subs] = odeToVectorField([DE1; DE2; DE3])
odesfcn = matlabFunction(VF, 'Vars',{'t','Y', 'c_1', 'c_2', 'k_1', 'k_2', 'k_n', 'k_nl', 'J_1', 'J_2', 'J_n'})
I leave the rest of the coding to you.
First, check to be certain all the parameters have been accounted for and declared in the syms call before you finish coding the differential equations symbolically.
Then use ‘odesfcn’ with ode45 or other appropriate numerical integration function of your choice as:
@(t,Y)odesfcn(t, Y, c_1, c_2, k_1, k_2, k_n, k_nl, J_1, J_2, J_n);
having assigned all the parameters in your code earlier.
.
  6 Commenti
Sam Butler
Sam Butler il 8 Feb 2023
@Walter Roberson and @Star Strider, thank you for your help with this. It looks like it was just the last few lines of code I was missing out on.
As I now have six equations, how can I numerically solve these? I would ideally like to plot some curves to show the system behaviour, for example torque v deflection and omega v time.
Thank you for your help!
Star Strider
Star Strider il 8 Feb 2023
@Walter Roberson is the hero here. I got reasonably far, however not far enough.
It would likely be best to use ode15s, since the parameter values span a wide range of orders-of-magnitude, and the system is likely ‘stiff’ as the result. Create a vector of appropriate initial conditions for the variables in ‘NewVar’ and then approach it as:
odesfcn = @(t,in2)[in2(5,:);in2(6,:);in2(7,:);in2(5,:).*(-9.372351067469629e+1)+in2(6,:).*9.098146252701012e+1+in2(7,:).*3.742048147686167-in2(1,:).*3.810889617686484e+6+in2(2,:).*3.810045490731264e+6+in2(3,:).*2.644126955220157e+3-in2(4,:).*1.8e+3-sqrt(-in2(1,:)+in2(3,:)).*2.591333333333333e+5+sqrt(-in2(5,:)+in2(7,:)).*9.69689410003742e+6;0.0;in2(5,:).*6.823609689525759e+1-in2(6,:).*6.823609689525759e+1+in2(1,:).*2.900034118048448e+6-in2(2,:).*2.900034118048448e+6;in2(5,:).*3.742048147686167e+1-in2(7,:).*3.742048147686167e+1+in2(1,:).*2.716726955220157e+4-in2(3,:).*2.716726955220157e+4-sqrt(-in2(5,:)+in2(7,:)).*9.69689410003742e+7] ;
ic = [...]; % Initial Conditions Vector
tspan = [0 5]; % Integration Time Limits
[t,y] = ode15s(odesfcn, ic, tspan);
figure
plot(t, y)
grid
legend(string(NewVar))
See the documentation on ode15s for details.
If you want the results reported at specific times, use:
tspan = linspace(0, 5, 100);
or something similar. The integration function will interpolate the results to those values.
.

Accedi per commentare.


Sam Chak
Sam Chak il 8 Feb 2023
Some of the the variable names in your original code are different from the ODEs shown in the image. So I fixed them and assigned some initial values. Also, I think that 800 rpm needs to be converted to rad/s.
tspan = linspace(0, 0.1, 10001);
x0 = [1 0.5 0.25 0 0 0 85.8702];
[t, x] = ode45(@odefcn, tspan, x0);
plot(t, x(:,1:3)), grid on, xlabel('t'), ylabel('\bf{\theta}')
legend('\theta_1', '\theta_2', '\theta_n')
function xdot = odefcn(t, x)
xdot = zeros(7, 1);
% Parameters
J_1 = 2.08e-4;
J_2 = 2.931e-3;
J_n = 4.0085e-4;
c_1 = 0.15;
c_2 = 0.2;
c_n = 0.015;
k_1 = 270;
k_2 = 8500;
k_n = 10.89;
k_nl = 38870;
% Motor Input Velocity
% Motor oscillates very fast, so T = linspace(0,100,5) is insufficient
d1theta_0 = 800*2*pi/60; % 800 rpm is not the same as 800 rad/s
d1theta_3 = 0.02*d1theta_0;
d1theta_4 = 0.005*d1theta_0;
f_mech = d1theta_0/(2*pi);
d1theta_out = d1theta_0 + d1theta_3*cos(72*pi*f_mech*t) + d1theta_4*cos(144*pi*f_mech*t);
% ODEs
xdot(1) = x(4); % x(1) is θ1, x(4) is θ1'
xdot(2) = x(5); % x(2) is θ2, x(5) is θ2'
xdot(3) = x(6); % x(3) is θn, x(6) is θn'
xdot(4) = (- (c_1*(x(4) - d1theta_out) + k_1*(x(1) - x(7)) - c_2*(x(5) - x(4)) - k_2*(x(2) - x(1)) - c_n*(x(6) - x(4)) - k_n*(x(3) - x(1)) - k_nl*(x(3) - x(1))^5))/J_1;
xdot(5) = (- (c_2*(x(5) - x(4)) + k_2*(x(2) - x(1))))/J_2;
xdot(6) = (- (c_n*(x(6) - x(4)) + k_n*(x(3) - x(1)) + k_nl*(x(3) - x(1))^5))/J_n;
xdot(7) = d1theta_out; % x(7) is θout
end

Tag

Prodotti


Release

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by