Solve a system of 3 second order equations

Question

Sam Butler 2023-2-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1908130-solve-a-system-of-3-second-order-equations

评论： Star Strider 2023-2-8

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Alan Stevens 2023-2-7

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1908130-solve-a-system-of-3-second-order-equations#answer_1165740

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

Star Strider 2023-2-7

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1908130-solve-a-system-of-3-second-order-equations#answer_1165750

在 MATLAB Online 中打开

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 个评论
显示 4更早的评论隐藏 4更早的评论

Sam Butler 2023-2-7

在 MATLAB Online 中打开

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 2023-2-7

在 MATLAB Online 中打开

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 2023-2-7

在 MATLAB Online 中打开

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]

Sam Butler 2023-2-8

@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 2023-2-8

@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 2023-2-8

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1908130-solve-a-system-of-3-second-order-equations#answer_1166430

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Solve a system of 3 second order equations

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（3 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

6 个评论
显示 4更早的评论隐藏 4更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Solve a system of 3 second order equations

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（3 个）

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

6 个评论 显示 4更早的评论隐藏 4更早的评论

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

6 个评论
显示 4更早的评论隐藏 4更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论