second order non-linear ODE and curve fitting: problem with initial conditions
2 views (last 30 days)
Show older comments
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');
0 Comments
Answers (1)
Torsten
on 8 Apr 2015
For theta_0=0,phi_dot_dummy=Inf for your initial conditions since you divide by sin^2(theta_0)=0.
Best wishes
Torsten.
2 Comments
Torsten
on 8 Apr 2015
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.
See Also
Categories
Find more on Assembly in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!