Re: Error Tracing - NDSolve
- To: mathgroup at smc.vnet.net
- Subject: [mg4122] Re: Error Tracing - NDSolve
- From: withoff (David Withoff)
- Date: Wed, 5 Jun 1996 01:39:03 -0400
- Organization: Wolfram Research, Inc.
- Sender: owner-wri-mathgroup at wolfram.com
In article <4obg9e$na9 at dragonfly.wolfram.com> Andrei Constantinescu
<constant at athena.polytechnique.fr> writes:
>
> Hi,
>
> is there any possibility in tracing an error in
> NDSolve ?
>
> I get the following message:
>
> In[2]:= sol = NDSolve[ Union[equation, condini] , {y1,y2,y3} , {t,
0., 20.},
> MaxSteps -> 2000 ]
>
> 1
> Power::infy: Infinite expression -- encountered.
> 0.
>
> ..... and unfortunately, the system seems to be well defined.
>
> Thanks,
>
> a + andrei
>
> PS I can send the entire file if there is a more precise interest
> in the question.
>
>
> ______________________________________________________________________
> Andrei Constantinescu constant at athena.polytechnique.fr
>
> LMS Ecole Polytechnique tel: (33)-1-69.33.33.30
> 91128 PALAISEAU cedex - FRANCE fax: (33)-1-69.33.30.26
>
>
Most difficulties with NDSolve fall into three categories:
1) Errors in the input, such as syntax errors, undefined symbols,
or types of equations (such as boundary value problems) that NDSolve
isn't (at least for now) designed to handle. I don't think the
problem you describe falls into this category.
2) Problems with preprocessing of the equations, such as in solving
for the initial conditions or the highest-order derivatives. I suspect
that your example falls into this category.
3) Difficulties that come up when the differential equation is solved,
such as singular, oscillatory, or unstable equations. If your claim
about the system being well-defined is correct, then your example
probably doesn't fall into this category.
Difficulties of type (1) are resolved by entirely mundane
procedures such as reading the manual, and since your example
doesn't appear to fall into that category I won't get into that.
There are several useful tricks for tracing through the rest of
the work that is done by NDSolve. One of the more useful tricks is
to define a function that always evaluates to zero, but that prints
or records useful information as a side effect, and to add this
function to one of your equations.
In[1]:= showprogress[p_?NumberQ] := (Print[p]; 0)
In[2]:= NDSolve[{f'[x] == f[x] + showprogress[x], f[1] == 1},
f[x], {x, 1, 2}]
1.
1.00141
1.00141
1.00283
1.00283
1.02341
1.02341
1.04399
1.04399
1.06457
1.06457
1.1238
1.1238
1.18304
1.18304
1.24228
1.24228
1.30152
1.30152
1.40002
1.40002
1.49853
1.49853
1.59703
1.59703
1.69553
1.69553
1.79403
1.79403
1.94641
1.94641
2.
2.
Out[2]= {{f[x] -> InterpolatingFunction[{1., 2.}, <>][x]}}
If you are impatient (like me) and your NDSolve example is big
and complicated, this approach gives you something to watch
while the calculation proceeds, and gives you reassurance that
NDSolve is making progress towards a solution.
In many cases this trick will also help you separate problems of
type (2) from problems of type (3).
In your example, if you see the mysterious Power::infy message before
you see the numbers printed by the showprogress function, then you
know that the problem must be in preprocessing of the equations,
which is a problem of type (2), and probably doesn't have anything
to do with the heart of the numerical differential equations solver,
or the fact that the equations are differential equations. If the
Power::infy appears after NDSolve starts solving the equations, then
the problem is probably somewhere in calculation of the derivatives,
or in something else related to construction of the solution.
Most problems of type (2) are in Solve, because the most complicated
part of preprocessing is done by Solve. NDSolve often asks Solve
to do exotic things that are perfectly legitimate, but that you would
probably never enter yourself.
You can find out what Solve was asked to do by adding a rule to Solve
that records its arguments as a side effect.
In[3]:= Unprotect[Solve]
Out[3]= {Solve}
In[4]:= Solve[p___] := $Failed /;
(Print["using Solve rule"]; AppendTo[args, {p}])
In[5]:= args = {}; NDSolve[{f'[x] == f[x], f[1] == 1}, f[x], {x, 1, 2}]
using Solve rule
using Solve rule
Out[5]= {{f[x] -> InterpolatingFunction[{1., 2.}, <>][x]}}
called once to solve for the initial conditions
In[6]:= args[[1]]
Out[6]= {{-1. + f[1.] == 0.}, {f[1.]}}
and once to solve for the highest-order derivatives
In[7]:= args[[2]]
Out[7]= {{-1. f[x] + f'[x] == 0}, {f'[x]}}
You can then apply Solve directly to these arguments, and thereby
check the most complicated step of preprocessing in isolation from
the rest of NDSolve. In the present example, I suspect that this
will not only expose the problem, but will also help you figure out
what, if anything, should be done about it. If my guess is correct,
the Power::infy message is harmless and can be ignored.
If this analysis leads to other questions, or you'd like someone to
take a look at your equations, you are of course welcome to send
them along.
Dave Withoff
Research and Development
Wolfram Research
==== [MESSAGE SEPARATOR] ====