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