Re: Limitations of NDSolve

*To*: mathgroup at smc.vnet.net*Subject*: [mg2540] Re: Limitations of NDSolve*From*: saarinen (Sirpa Saarinen)*Date*: Fri, 17 Nov 1995 00:20:11 -0500*Organization*: Wolfram Research, Inc.

In article <48c046$3f2 at ralph.vnet.net> David Garza <dgarza at utsi.edu> writes: > 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 > > > I tried this problem with V2.1 and V2.2 and I cannot re-produce the problem that you describe. Maybe you could tell us which version you use, which environment you have and what else is different from the system below (with V2.2). > > In[1]:= 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 work s *) > al=0; (* Angle of attack, set to zero to see if this works *) > rho=.00276; (* Air density, try a constant for now *) > > > > In[2]:= > In[3]:= > In[4]:= > In[5]:= > In[6]:= > In[7]:= > In[8]:= > In[9]:= > In[10]:= > In[11]:= > In[12]:= > In[13]:= In[13]:= In[13]:= NDSolve[{ > 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]), > In[4]:= > In[5]:= > In[6]:= > In[7]:= > In[8]:= > In[9]:= > In[10]:= > In[11]:= > In[12]:= > In[13]:= In[13]:= In[13]:= NDSolve[{ > 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) ,r[0]==r0, > th[0]==.1, > v[0]==.1, > gam[0]==Pi*89.995/180, > m[0]==10000 > } > ,{r,th,v,gam,m},{t ,0,100}] > > > Out[13]= {{r -> InterpolatingFunction[{0., 100.}, <>], > > > th -> InterpolatingFunction[{0., 100.}, <>], > > > v -> InterpolatingFunction[{0., 100.}, <>], > > > gam -> InterpolatingFunction[{0., 100.}, <>], > > > m -> InterpolatingFunction[{0., 100.}, <>]}} > > In[14]:= In[14]:= > > > Sirpa Saarinen > saarinen at wri.com