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] ====