Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
1995
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 1995

[Date Index] [Thread Index] [Author Index]

Search the Archive

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



  • Prev by Date: Re: turn off symbols in MultipleListPlot?
  • Next by Date: Why no simplification ?
  • Previous by thread: Limitations of NDSolve
  • Next by thread: Final comments on CoefficientList