second order non-linear ODE and curve fitting: problem with initial conditions
1 次查看(过去 30 天)
显示 更早的评论
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 个评论
回答(1 个)
Torsten
2015-4-8
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 个评论
Torsten
2015-4-8
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.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Assembly 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!