Solve a system of 3 second order equations

조회 수: 1 (최근 30일)
Sam Butler
Sam Butler 2023년 2월 7일
댓글: 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.

답변 (3개)

Alan Stevens
Alan Stevens 2023년 2월 7일
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 2023년 2월 7일
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
Sam Butler
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
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.
.

댓글을 달려면 로그인하십시오.


Sam Chak
Sam Chak 2023년 2월 8일
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

카테고리

Help CenterFile Exchange에서 Programming에 대해 자세히 알아보기

태그

제품


릴리스

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by