       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==r0,
>   th==.1,
>   v==.1,
>   gam==Pi*89.995/180,
>   m==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:= 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:=
>      In:=
>      In:=
>      In:=
>      In:=
>      In:=
>      In:=
>      In:=
>      In:=
>      In:=
>      In:=
>      In:= In:= In:= 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:=
>      In:=
>      In:=
>      In:=
>      In:=
>      In:=
>      In:=
>      In:=
>      In:=
>      In:= In:= In:= 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==r0,
>        th==.1,
>        v==.1,
>        gam==Pi*89.995/180,
>        m==10000
>      }
>
,{r,th,v,gam,m},{t
,0,100}]
>
>
>      Out= {{r -> InterpolatingFunction[{0., 100.}, <>],
>
>      >     th -> InterpolatingFunction[{0., 100.}, <>],
>
>      >     v -> InterpolatingFunction[{0., 100.}, <>],
>
>      >     gam -> InterpolatingFunction[{0., 100.}, <>],
>
>      >     m -> InterpolatingFunction[{0., 100.}, <>]}}
>
>      In:= In:=
>
>
>      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