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

Question

Martin Cheung 2015-4-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/196700-second-order-non-linear-ode-and-curve-fitting-problem-with-initial-conditions

评论： Torsten 2015-4-8

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