Solve a system of 3 second order equations

1 view (last 30 days)
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.

Answers (3)

Alan Stevens
Alan Stevens on 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 on 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.
Star Strider
Star Strider on 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);
plot(t, y)
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.

Sign in to comment.

Sam Chak
Sam Chak on 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





Community Treasure Hunt

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

Start Hunting!