![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1287650/image.png)
Solve a system of 3 second order equations
2 次查看(过去 30 天)
显示 更早的评论
I want to solve the following 3 equations using MATLAB:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1287540/image.png)
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 个评论
回答(3 个)
Alan Stevens
2023-2-7
Write your equations as 6 first order equations like so:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1287650/image.png)
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 个评论
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 个评论
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))
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
2023-2-8
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 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!