second order non-linear ODE and curve fitting: problem with initial conditions

I am trying to solve a 2nd order nonlinear ODE. This problem arises from the rotation of a symmetric rigid body (e.g. spinning a top) with 1 endpoint pinned to ground.

Theta = angle between the line to centre of mass of the body and the usual z-axis (the inclination of the body).

Phi = rotation of body about z-axis.

Psi = rotaton of the body about its symmetry axis (the spin of the body itself).

I need to sketch Theta(t) vs t for a given set of initial conditions.

This code works fine as long as Theta_0 ~= 0.

When Theta_0 == 0, y1 becomes NaN and the program cannot proceed to the next section (curve fitting).

I cannont see why this is happening.

Is there any way to solve this problem?

if true
  % code
end
%
%
%
% Set initial conditions.
theta_0 = input('Enter intial value of Theta (from 0 to pi) : ');
theta_dot_0 = input('Enter intial value Theta_dot: ');
phi_0 = input('Enter intial value of Phi (from 0 to 2*pi) : '); 
phi_dot_0 = input('Enter intial value of Phi_dot: '); 
psi_0 = input('Enter initial value of Psi (from 0 to 2*pi) : '); 
psi_dot_0 = input('Enter initial value of Psi_dot: ');
C = input('Enter C (C>0) : ');
%
%
%
alpha = C*(psi_dot_0 + phi_dot_0*cos(theta_0)); 
beta = alpha*cos(theta_0) + phi_dot_0*sin(theta_0)^2;
%
%
%
syms theta(t)
phi_dot_dummy = (beta - alpha*cos(theta))/sin(theta)^2; 
V1 = odeToVectorField(diff(theta,2) == (phi_dot_dummy*(phi_dot_dummy*cos(theta) - alpha) + 1)*sin(theta));
M1 = matlabFunction(V1,'vars', {'t','Y'}); 
[x1,y1] = ode45(M1,[0 20],[theta_0 theta_dot_0]);
%
%
%
theta_fit = fit(x1,y1(:,1),'fourier2');

Answer 1

[Torsten]

For theta_0=0,phi_dot_dummy=Inf for your initial conditions since you divide by sin^2(theta_0)=0.

Best wishes

Torsten.

[Martin Cheung]

Thanks!

I figured if I set theta_0 = 1e-20 when the input of theta_0 = 0, the model seems to behave correctly.

But is this the correct way to do it? By making theta_0 small?

[Torsten]

If you calculate phi_dot_dummy at t=0, you'll see that the actual result is phi_dot_0. Thus you can circumvent y1=NaN by setting

phi_dot_dummy=phi_dot_0 for t=0 and

phi_dot_dummy=(beta - alpha*cos(theta))/sin(theta)^2 else.

Best wishes

Torsten.