Limitations of NDSolve
- Subject: [mg2515] Limitations of NDSolve
- From: dgarza at utsi.edu (David Garza)
- Date: Wed, 15 Nov 1995 07:01:42 GMT
- Approved: usenet@wri.com
- Distribution: local
- Newsgroups: wri.mathgroup
- Organization: University of Tennessee, Knoxville
- Sender: daemon at wri.com ( )
I have tried using Mathematica to numerically solve the equations of
motion for a vehicle in flight about a round Earth. The flight is in
two dimensions and the equations consist of a set of five nonlinear
ordinary differential equations. I have tried the same thing without
problems using computer programs and Excel. Even the Runge-Kutta
package supplied with Mathematica does fine. Unfortunately, NDSolve
will not work, instead reporting an infinite expression. The initial
conditions normally provide a small, nonzero value for the velocity
to eliminate the problem of velocity in the denominator in the fourth
equation. This has been something I've tried off and on a couple of
times in the past, and in talking with others I have found that NDSolve
has had trouble with other equations when those equations have not
presented problems to other means of computing numerical solutions,
like writing a program in Fortran or C. What sort of limitations are
there on NDSolve with regard to the problems it can handle? The system
of equations I use are included below, along with initial conditions and
the parameters necessary for solution.
ODE's:
{
th'[t]==-(v[t] Cos[gam[t]]/r[t]),
r'[t]==v[t] Sin[gam[t]],
v'[t]==(-drag+thrust Cos[al+eps]-g m[t]*
Sin[gam[t]])/m[t],
gam'[t]==(lift-m[t] g Cos[gam[t]]+ m[t]*
v[t]^2 Cos[gam[t]]/r[t]+thrust Sin[al+eps])/
(m[t]*v[t]),
m'[t]==-thrust/(Isp g)
}
Initial Conditions:
{
r[0]==r0,
th[0]==.1,
v[0]==.1,
gam[0]==Pi*89.995/180,
m[0]==10000
}
Other important functions:
thrust=1.3*m0*g;
drag=.5*rho*v[t]^2*area*cd;
lift=0; (* This is a rocket, so there is little lift - neglect it now *)
g=9.81; (* surface gravitational acceleration *)
m0=10000; (* Initial mass *)
r0=6376000; (* Radius of the Earth *)
area=10; (* Rocket frontal area *)
cd=.1; (* Drag coefficient *)
Isp=400; (* Engine specific impulse *)
eps=0; (* Engine gimbal angle, try 0 for now to see if this works *)
al=0; (* Angle of attack, set to zero to see if this works *)
rho=.00276; (* Air density, try a constant for now *)
Thanks,
--
David Garza
dgarza at utsi.edu