Re: solving differential equations
- To: mathgroup at smc.vnet.net
- Subject: [mg3544] Re: solving differential equations
- From: Robert Knapp <rknapp>
- Date: Thu, 21 Mar 1996 01:10:25 -0500
- Organization: Wolfram Research, Inc.
- Sender: owner-wri-mathgroup at wolfram.com
Steven Wilkinson wrote: > > Mathgroup, > > I have a system of ODE I wish to solve numerically of the form X'[t] = > F(X[t])/|F(X[t])| for a given function > F : plane -> plane > and given intial conditions. NDSolve works well on this except where the > system is ill-behaved, in particular, where F = 0. I can replace this > system with another one that has equivalent solutions and is usually > well-behaved where F = 0. However, I don't want to make the replacement > unless I have to, and only near where F = 0. (I wish to construct a command > that gives the solution for an inputted function F.) > > Is there any way possible to do this with NDSolve? In other words, as it > steps along have it check the value of F. If its magnitude ever gets too > close to 0, have it use the other system. > > If that is not possible, does anyone have a solver package that I could > modify to do the above? > > Thanks, > Steve Wilkinson > Dept of Math & CS > Northern Kentucky University > wilkinson at nku.edu The easiest thing to do is to redefine the function so that you change to the equivalent solution only when F is near zero. Below is a very simple example of how you can do this. In[1]:= NDSolve[{x'[t] == t/Sin[t],x[0] == -1},x,{t,0,1}] Infinity::indet: Indeterminate expression 0. 1. ComplexInfinity encountered. NDSolve::ndnum: Differential equation does not evaluate to a number at t = 0.. Out[1]= {{x -> InterpolatingFunction[{0., 0.}, <>]}} In[2]:= Series[t/Sin[t],{t,0,3}] Out[2]= (* find an approximation near the removeable singularity *) 2 t 4 1 + -- + O[t] 6 In[3]:= (* define the function, so that the approximation is used near the removable singularity. 100 $MachineEpsilon was an arbitrary choice, but I checked it with a plot to see that it matched the approximation well *) f[t_] := If[Abs[Sin[t]] < 100 $MachineEpsilon, 1+t^2/6,t/Sin[t]] In[4]:= NDSolve[{x'[t] == f[t],x[0] == 1},x,{t,-1,1}] Out[4]= {{x -> InterpolatingFunction[{-1., 1.}, <>]}} Note that you case is slightly different since your right hand side has the function with singularity depending on the dependent variable. You can still use a similar definition and use f[x[t]] instead of just f[t]. Rob Knapp Wolfram Research, Inc. http://www.wri.com/~rknapp ==== [MESSAGE SEPARATOR] ====