Solving system of nonlinear differential equations using ode45

3 vues (au cours des 30 derniers jours)
Miroslav Jovic
Miroslav Jovic le 11 Juil 2024
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 commentaires
Torsten
Torsten le 11 Juil 2024
sin(0) = 0 - thus you divide by 0.
Miroslav Jovic
Miroslav Jovic le 17 Juil 2024
You are right, that is one of the problems. Thanks!

Connectez-vous pour commenter.

Réponses (3)

Francisco J. Triveno Vargas
Francisco J. Triveno Vargas le 12 Juil 2024
Modifié(e) : Torsten le 12 Juil 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)
  1 commentaire
Miroslav Jovic
Miroslav Jovic le 17 Juil 2024
Looks like initial values are another problem, thank you!

Connectez-vous pour commenter.


Sam Chak
Sam Chak le 11 Juil 2024
Modifié(e) : Sam Chak le 11 Juil 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 commentaires
Miroslav Jovic
Miroslav Jovic le 17 Juil 2024
I have tried for constant values at first, and now i've calculate real torques and put it in code. Thank you for your help!
Sam Chak
Sam Chak le 17 Juil 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!

Connectez-vous pour commenter.


Muhammed abdulmalek
Muhammed abdulmalek le 11 Juil 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')

Catégories

En savoir plus sur Particle & Nuclear Physics dans Help Center et File Exchange

Produits


Version

R2015a

Community Treasure Hunt

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

Start Hunting!

Translated by