Solving system of nonlinear differential equations using ode45

1 次查看（过去 30 天）
Miroslav Jovic
Miroslav Jovic 2024-7-11
编辑： Walter Roberson 2024-7-17
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.

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

回答（3 个）

Francisco J. Triveno Vargas
编辑：Torsten 2024-7-12
Ran in:
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
Sam Chak 2024-7-11
编辑：Sam Chak 2024-7-11
Ran in:
Hi @Miroslav Jovic
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 个评论
Miroslav Jovic
Miroslav Jovic 2024-7-17
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 2024-7-17
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
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')

类别

在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

标签

产品


版本

R2015a

Community Treasure Hunt

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

Start Hunting!

Translated by