Solve a system of 3 second order equations

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.

Answer 1

[Alan Stevens]

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.

Answer 2

[Star Strider]

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.

.

[Sam Butler]

Thank you for your help.

The input velocity (from a motor), is characterised by the following equation

in the interval T = 0 to 5seconds

d1theta_0 can be inputted on each run, e.g. 800rpm.

I have attempted to do this in the code however I also need the integral for the theta_out term.

Could you recomend a way of doing this please? I have attached the completed code so far.

I am new to MATLAB so appreciate nay help that I can get.

syms c_1 c_2 c_n 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 

%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;

%Define Variables
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);

%Motor Input Vel
d1theta_0 = 800;
d1theta_3 = 0.02*d1theta_0;
d1theta_4 = 0.005*d1theta_0;
f_mech = d1theta_0/(2*pi);
T = linspace(0,100,5);
d1theta_out = d1theta_0 + d1theta_3*cos(72*pi*f_mech*T) + d1theta_4*cos(144*pi*f_mech*T);
% Add theta_out equation here.

%Equations of Motion
DE1 = J_1*d2theta_1 + c_1*(d1theta1 - d1theta_out) + k_1*(theta_1 - theta_out)...
    - c_2*(d2theta_2 - d2theta_1) - k_2*(theta_2 - theta_1)...
    - c_n*(d2theta_n - d2theta_1) - k_n*(theta_n - theta_1)...
    - k_nl*(theta_n - theta_1).^5 = 0;
DE2 = J_n*d2theta_n + c_n*(d1theta_n - d1theta_1) + k_n*(theta_n - theta_1)...
    + k_nl*(d1theta_n - d1theta_1).^5 = 0;
DE3 = J_2*d2theta_2 + c_2*(d1theta_2 - d1theta_1) + k_2*(theta_2 - theta_1) = 0


[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'})

[Star Strider]

I tried several different approaches (odeToVectorField, reduceDifferentialOrder, odeFunction, and others), and I cannot get it into a form that the numerical integrators will accept.

syms c_1 c_2 c_n 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 
sympref('AbbreviateOutput',false);
%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;

%Define Variables
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);

%Motor Input Vel
d1theta_0 = 800;
d1theta_3 = 0.02*d1theta_0;
d1theta_4 = 0.005*d1theta_0;
f_mech = d1theta_0/(2*pi);
T = linspace(0,100,5);
d1theta_out = d1theta_0 + d1theta_3*cos(72*pi*f_mech*T) + d1theta_4*cos(144*pi*f_mech*T);
% Add theta_out equation here.

%Equations of Motion
DE1 = J_1*d2theta_1 == -(c_1*(d1theta_1 - d1theta_out) + k_1*(theta_1 - theta_out)...
    - c_2*(d2theta_2 - d2theta_1) - k_2*(theta_2 - theta_1)...
    - c_n*(d2theta_n - d2theta_1) - k_n*(theta_n - theta_1)...
    - k_nl*(theta_n - theta_1).^5);
DE2 = J_n*d2theta_n == -(c_n*(d1theta_n - d1theta_1) + k_n*(theta_n - theta_1)...
    + k_nl*(d1theta_n - d1theta_1).^5);
DE3 = J_2*d2theta_2 == -(c_2*(d1theta_2 - d1theta_1) + k_2*(theta_2 - theta_1));


[NewEQ,NewVar] = reduceDifferentialOrder([DE1,DE2,DE3], [theta_1, theta_2,theta_n,theta_out])

NewEQ = 

NewVar = 

% odesfcn = odeFunction(NewEQ, NewVar)

% [VF, Subs] = odeToVectorField(DE1, DE2, DE3)
odesfcn = matlabFunction(NewEQ)%, 'Vars',{'t','Y', 'c_1', 'c_2', 'k_1', 'k_2', 'k_n', 'k_nl', 'J_1', 'J_2', 'J_n'})

Warning: Function 'Dtheta_1t' not verified to be a valid MATLAB function.

Warning: Function 'theta_1' not verified to be a valid MATLAB function.

Warning: Function 'theta_2' not verified to be a valid MATLAB function.

Warning: Function 'theta_n' not verified to be a valid MATLAB function.

Warning: Function 'theta_out' not verified to be a valid MATLAB function.

Warning: Function 'Dtheta_2t' not verified to be a valid MATLAB function.

Warning: Function 'Dtheta_nt' not verified to be a valid MATLAB function.

odesfcn = function_handle with value:
    @(t)[Dtheta_1t(t).*(3.0./2.0e+1)+theta_1(t).*8.78089e+3-theta_2(t).*8.5e+3-theta_n(t).*1.089e+1-theta_out(t).*2.7e+2+diff(Dtheta_1t(t),t).*2.15208e-1-diff(Dtheta_2t(t),t)./5.0-diff(Dtheta_nt(t),t).*(3.0./2.0e+2)+(theta_1(t)-theta_n(t)).^5.*3.887e+4-1.23e+2;Dtheta_1t(t).*(3.0./2.0e+1)+theta_1(t).*8.78089e+3-theta_2(t).*8.5e+3-theta_n(t).*1.089e+1-theta_out(t).*2.7e+2+diff(Dtheta_1t(t),t).*2.15208e-1-diff(Dtheta_2t(t),t)./5.0-diff(Dtheta_nt(t),t).*(3.0./2.0e+2)+(theta_1(t)-theta_n(t)).^5.*3.887e+4-1.182055469163792e+2;Dtheta_1t(t).*(3.0./2.0e+1)+theta_1(t).*8.78089e+3-theta_2(t).*8.5e+3-theta_n(t).*1.089e+1-theta_out(t).*2.7e+2+diff(Dtheta_1t(t),t).*2.15208e-1-diff(Dtheta_2t(t),t)./5.0-diff(Dtheta_nt(t),t).*(3.0./2.0e+2)+(theta_1(t)-theta_n(t)).^5.*3.887e+4-1.218218161547023e+2;Dtheta_1t(t).*(3.0./2.0e+1)+theta_1(t).*8.78089e+3-theta_2(t).*8.5e+3-theta_n(t).*1.089e+1-theta_out(t).*2.7e+2+diff(Dtheta_1t(t),t).*2.15208e-1-diff(Dtheta_2t(t),t)./5.0-diff(Dtheta_nt(t),t).*(3.0./2.0e+2)+(theta_1(t)-theta_n(t)).^5.*3.887e+4-1.185732137846494e+2;Dtheta_1t(t).*(3.0./2.0e+1)+theta_1(t).*8.78089e+3-theta_2(t).*8.5e+3-theta_n(t).*1.089e+1-theta_out(t).*2.7e+2+diff(Dtheta_1t(t),t).*2.15208e-1-diff(Dtheta_2t(t),t)./5.0-diff(Dtheta_nt(t),t).*(3.0./2.0e+2)+(theta_1(t)-theta_n(t)).^5.*3.887e+4-1.196183945125302e+2;Dtheta_1t(t).*(-3.0./2.0e+2)+Dtheta_nt(t).*(3.0./2.0e+2)-theta_1(t).*1.089e+1+theta_n(t).*1.089e+1-(Dtheta_1t(t)-Dtheta_nt(t)).^5.*3.887e+4+diff(Dtheta_nt(t),t).*4.0085e-4;Dtheta_1t(t).*(-1.0./5.0)+Dtheta_2t(t)./5.0-theta_1(t).*8.5e+3+theta_2(t).*8.5e+3+diff(Dtheta_2t(t),t).*2.931e-3;Dtheta_1t(t)-diff(theta_1(t),t);Dtheta_2t(t)-diff(theta_2(t),t);Dtheta_nt(t)-diff(theta_n(t),t)]

% odesfcn = matlabFunction(NewEQ)

It may not be possible to solve this.

.

[Walter Roberson]

I was getting misled by the way the variable names were appearing in the equations, so I changed variable names to ones that MATLAB would know how to render appropriate.

I also followed the massMatrix logic in the odeFunction example, and that seemed to have worked.

syms c_1 c_2 c_n 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 
sympref('AbbreviateOutput',false);
Q = @(v) sym(v);
%Parameters
J_1 = Q(208)*sym(10)^(-6);
J_2 = Q(2931)*sym(10)^(-6);
J_n = Q(40085)*sym(10)^(-8)

J_n = 

c_1 = Q(15)*sym(10)^(-2);
c_2 = Q(2)*sym(10)^(-1);
c_n = Q(15)*sym(10)^(-3);
k_1 = Q(270);
k_2 = Q(8500);
k_n = Q(1089)*sym(10)^(-2);
k_nl = Q(38870);
Pi = sym(pi);

%Define Variables
theta_1_prime = diff(theta_1);
theta_1_dprime = diff(theta_1_prime);
theta_2_prime = diff(theta_2);
theta_2_dprime = diff(theta_2_prime);
theta_n_prime = diff(theta_n);
theta_n_dprime = diff(theta_n_prime);
theta_out_prime = diff(theta_out);

%Motor Input Vel
theta_0_prime = Q(800);
theta_3_prime = 0.02*theta_0_prime;
theta_4_prime = 0.005*theta_0_prime;
f_mech = theta_0_prime/(2*Pi);
T = linspace(0,100,5);
d1theta_out = theta_0_prime + theta_3_prime*cos(72*Pi*f_mech*T) + theta_4_prime*cos(144*Pi*f_mech*T);
% Add theta_out equation here.

%Equations of Motion
DE1 = J_1*theta_1_dprime == -(c_1*(theta_1_prime - theta_out_prime) + k_1*(theta_1 - theta_out)...
    - c_2*(theta_2_dprime - theta_1_dprime) - k_2*(theta_2 - theta_1)...
    - c_n*(theta_n_dprime - theta_1_dprime) - k_n*(theta_n - theta_1)...
    - k_nl*sqrt(theta_n - theta_1));
DE2 = J_n*theta_n_dprime == -(c_n*(theta_n_prime - theta_1_prime) + k_n*(theta_n - theta_1)...
    + k_nl*sqrt(theta_n_prime - theta_1_prime));
DE3 = J_2*theta_2_dprime == -(c_2*(theta_2_prime - theta_1_prime) + k_2*(theta_2 - theta_1));
eqns = [DE1, DE2, DE3].'

eqns(t) = 

vars = findSymType(eqns, 'symfun');

[NewEQ,NewVar] = reduceDifferentialOrder(eqns, vars)

NewEQ = 

NewVar = 

[M,F] = massMatrixForm(NewEQ,NewVar)

M = 

F = 

f = M\F

Warning: Solution is not unique because the system is rank-deficient.

f = 

odesfcn = odeFunction(f, NewVar)

odesfcn = function_handle with value:
    @(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]

[Star Strider]

@Walter Roberson — Thank you!

I figured that the system has a solution, however it was completely eluding me. I saw massMatrixForm in the other documentation, however it did not occur to me to apply it here.

[Sam Butler]

@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]

@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.

.

Answer 3

[Sam Chak]

Hi @Sam Butler

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