Solving system of nonlinear differential equations using ode45

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.

[Torsten]

sin(0) = 0 - thus you divide by 0.

[Miroslav Jovic]

You are right, that is one of the problems. Thanks!

Answer 1

[Francisco J. Triveno Vargas]

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)

[Miroslav Jovic]

Looks like initial values are another problem, thank you!

Answer 2

[Sam Chak]

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)

[Miroslav Jovic]

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]

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!

Answer 3

[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')