Solving system of nonlinear differential equations using ode45
2 views (last 30 days)
Show older comments
I have a question for system of ordinary differential equations, because Matlab gives some strange solution as output. There is the code:
Jo1=1;
Jo2=2;
Jo3=3;
Mo1=1;
Mo2=1;
Mo3=1;
f=@(t,x)[x(4).*sin(x(3))./sin(x(2))+x(5).*cos(x(3))./sin(x(2));
cos(x(3)).*x(4)-sin(x(3)).*x(5);
-x(4).*sin(x(3)).*cos(x(2))./sin(x(2))-x(6).*cos(x(3)).*cos(x(2))./sin(x(2));
(Mo1-(Jo3-Jo2).*x(6).*x(5))./Jo1;
(Mo2-(Jo1-Jo3).*x(4).*x(6))./Jo2;
(Mo3-(Jo2-Jo1).*x(5).*x(4))./Jo3];
[t, x]= ode45(f, [0,1],[0,0,0,0,0,0]);
I get solution for x4, x5 and x6, but for x1, x2 and x3 the solution is NaN.
So if someone have any advice it would help.
2 Comments
Answers (3)
Francisco J. Triveno Vargas
on 12 Jul 2024
Edited: Torsten
on 12 Jul 2024
A simples tests is change the initial values like this:
clear
close all
Jo1=1;
Jo2=2;
Jo3=3;
Mo1=1;
Mo2=1;
Mo3=1;
f=@(t,x)[x(4).*sin(x(3))./sin(x(2))+x(5).*cos(x(3))./sin(x(2));
cos(x(3)).*x(4)-sin(x(3)).*x(5);
-x(4).*sin(x(3)).*cos(x(2))./sin(x(2))-x(6).*cos(x(3)).*cos(x(2))./sin(x(2));
(Mo1-(Jo3-Jo2).*x(6).*x(5))./Jo1;
(Mo2-(Jo1-Jo3).*x(4).*x(6))./Jo2;
(Mo3-(Jo2-Jo1).*x(5).*x(4))./Jo3];
[t, x]= ode45(f, [0,1],[0.01,0.01,0.01,0,0,0]);
figure(10)
plot(t,x)
Sam Chak
on 11 Jul 2024
Edited: Sam Chak
on 11 Jul 2024
It appears that your current formulations may be erroneous. If you are not comfortable memorizing the required formulas, I would suggest referring to established textbooks to copy the most commonly used and accepted expressions.
Additionally, since you are applying constant torques of Mo1 = 1, Mo2 = 1, and Mo3 = 1 to the system, the system is not going to remain in the desired equilibrium point of [0, 0, 0, 0, 0, 0] (also the initial values). This may indicate the need to further review your system dynamics and control strategies.
Rotational kinematics in terms of a 3-2-1 Euler rotation sequence:
Here is a simple demo:
%% The System
function dx = ode(t, x)
% the parameters
Jo1 = 1;
Jo2 = 2;
Jo3 = 3;
% the inputs (need to be designed)
Mo1 = - 1*x(1) - 2*x(4);
Mo2 = - 1*x(2) - 2*x(5);
Mo3 = - 1*x(3) - 2*x(6);
% the differential equations
dx = [x(4) + sin(x(1))*tan(x(2))*x(5) + cos(x(1))*tan(x(2))*x(6)
cos(x(1))* x(5) - sin(x(1))* x(6)
sin(x(1))*sec(x(2))*x(5) + cos(x(1))*sec(x(2))*x(6)
(Mo1 + (Jo2 - Jo3)*x(2)*x(3))/Jo1
(Mo2 + (Jo3 - Jo1)*x(3)*x(1))/Jo2
(Mo3 + (Jo1 - Jo2)*x(1)*x(2))/Jo3];
end
%% Run the simulation
tspan = [0, 15];
x0 = deg2rad([3, 6, 9, 0, 0, 0]);
[t, x] = ode45(@ode, tspan, x0);
plot(t, rad2deg(x(:,1:3))), grid on, xlabel('t'), ylabel('Angle (degree)')
title('Time responses of Euler angles')
legend('\theta_{1}', '\theta_{2}', '\theta_{3}', 'fontsize', 14)
2 Comments
Sam Chak
on 17 Jul 2024
Hi @Miroslav Jovic, I'm glad to hear that. If you found the example and code helpful, please consider clicking 'Accept' ✔️ on the answer. Additionally, you can show your appreciation by voting 👍 other answers as a token of support for knowledge sharing. Your support is greatly appreciated!
Muhammed abdulmalek
on 11 Jul 2024
Jo1 = 1;
Jo2 = 2;
Jo3 = 3;
Mo1 = 1;
Mo2 = 1;
Mo3 = 1;
f = @(t,x) [x(4) + sin(x(1)*tan(x(2))*x(5) + cos(x(1))*tan(x(2))*x(6));
cos(x(1)*x(5) - sin(x(1))*x(6));
sin(x(1)*sn(x(2))*x(5) + cos(x(1))*sn(x(2))*x(6));
(Mo1 - (Jo3 - Jo2)*x(3)*x(2))/Jo1;
(Mo2 - (Jo1 - Jo3)*x(1)*x(3))/Jo2;
(Mo3 - (Jo2 - Jo1)*x(2)*x(1))/Jo3];
[t, x]= ode45(f, [0, 1], [0, 0, 0, 0, 0, 0]);
plot(t, x), ızgara açık, xlabel('t')
0 Comments
See Also
Categories
Find more on Particle & Nuclear Physics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!