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