NDSolve
Usage Message: NDSolve[eqns, y, {x, xmin, xmax}] finds a numerical solution to the ordinary differential equations
eqns for the function y with the independent variable x in the range xmin to xmax. NDSolve[eqns, y, {x, xmin, xmax}, {t, tmin, tmax}] finds a numerical
solution to the partial differential equations eqns. NDSolve[eqns, {y1, y2, ... }, {x, xmin, xmax}] finds numerical solutions for the functions yi.
Attributes[NDSolve] = {Protected} Options: AccuracyGoal
-> Automatic
Compiled
-> True
DifferenceOrder
-> Automatic
InterpolationPrecision
-> Automatic
MaxRelativeStepSize
-> 1
MaxSteps
-> Automatic
MaxStepSize
-> Infinity
Method
-> Automatic
PrecisionGoal
-> Automatic
SolveDelayed
-> False
StartingStepSize
-> Automatic
StoppingTest
-> None
WorkingPrecision
-> 16
Notes: InterpolatingFunction Expressions And Non-Numerical Coefficients There are
two known classes of examples in which NDSolve will incorrectly generate an error message indicating that the equations are not numerical. Both classes
of examples involve incorrect use of complex numbers within the NDSolve function. A simple workaround which solves the problem in most examples is to
use Re
to discard the imaginary part of the variable in all places in the input where the variable must be real. Although the most commonly reported examples of
these problems involve InterpolatingFunction
expressions, the underlying problems are not related to InterpolatingFunction
, and can come up when any function that does not accept complex numbers is used in the input to NDSolve. The first class of examples involves the use of explicit complex numbers in the input. Here is a typical example.
In[1]:= c = Interpolation[{1, 2, 3, 4, 5}]
Out[1]= InterpolatingFunction[{{1, 5}}, <>]
In[2]:= NDSolve[{f'[x] == I c[x], f[1] == 0}, f, {x, 1, 5}]
NDSolve::ndnum:
The right-hand side of the differential equation does not evaluate to a
number at x == 1..
Out[2]= NDSolve[{f'[x] == I InterpolatingFunction[{{1, 5}}, <>][x],
> f[1] == 0}, f, {x, 1, 5}]
The problem in the example above is related to the use of the Compile
function, which is called automatically from NDSolve. When the input includes complex numbers, the Compile
function uses complex values for all of the numerical expressions in the input, including the variable x. Although the imaginary part of x is
always zero, InterpolatingFunction
will not give a numerical result unless the argument is real. This behavior of InterpolatingFunction
can be seen by evaluating c[1.0 + 0.0 I] directly, outside of the NDSolve function.
In[3]:= c[1.0]
Out[3]= 1.
In[4]:= c[1.0 + 0.0 I]
Out[4]= InterpolatingFunction[{{1, 5}}, <>][1. + 0. I]
You can correct the problem in this example either by turning off the compiler using the Compiled
-> False
option in NDSolve, or by replacing c[x] with c[Re[x]]. The preferred approach in other examples will depend on the details of the example.
Out[4]= InterpolatingFunction[{{1, 5}}, <>][1. + 0. I]
In[5]:= NDSolve[{f'[x] == I c[x], f[1] == 0}, f, {x, 1, 5},
Compiled -> False]
Out[5]= {{f -> InterpolatingFunction[{{1., 5.}}, <>]}}
In[6]:= NDSolve[{f'[x] == I c[Re[x]], f[1] == 0}, f, {x, 1, 5}]
Out[6]= {{f -> InterpolatingFunction[{{1., 5.}}, <>]}}
The second class of examples in which NDSolve may report that the equations are not numerical involves the use of InterpolatingFunction
expressions in boundary value problems. In these examples, the unnecessary complex numbers are generated within the algorithm that is used to solve boundary value problems, rather than within the compiler. Here is a typical example.
In[7]:= d = Interpolation[{1, 1, 1, 1, 1}]
Out[7]= InterpolatingFunction[{{1, 5}}, <>]
In[8]:= NDSolve[{d[x] f''[x] == -1, f[1] == 0, f[5] == 0}, f, {x, 1, 5}]
NDSolve::pcnan:
Coefficients of the differential equation are not numbers or only one
linear nth order ordinary differential equation can be solved.
Out[8]= NDSolve[{InterpolatingFunction[{{1, 5}}, <>][x] f''[x] == -1,
> f[1] == 0, f[5] == 0}, f, {x, 1, 5}]
You can correct the problem in this example by replacing d[x] with d[Re[x]].
In[9]:= NDSolve[{d[Re[x]] f''[x] == -1, f[1] == 0, f[5] == 0}, f, {x, 1, 5}]
Out[9]= {{f -> InterpolatingFunction[{{1., 5.}}, <>]}}
Both of these problems are being investigated for the next version of Mathematica Inconsistent Boundary and Initial Conditions If the boundary and initial conditions for a partial differential equation are not consistent with each other, NDSolve will generate a warning message, and the calculation will ignore at least one of the conditions. You can solve this problem by changing your input so that the boundary and initial conditions are no longer inconsistent. Here is a typical example in which the initial condition is inconsistent with one of the boundary conditions. The initial condition y[x,0]==0 indicates that y[0,0] is 0, while the boundary condition y[0,t]==1 indicates that y[0,0] is 1. In this example, after generating a warning message, the calculation effectively uses y[0,t]==0 as the boundary condition.
In[1]:= NDSolve[{D[y[x, t], t] == D[y[x, t], x, x],
y[x, 0] == 0, y[0, t] == 1,
Derivative[1, 0][y][1, t] == 0},
y[x, t], {x, 0, 1}, {t, 0, 1}]
NDSolve::ibcinc: Warning: Boundary and initial conditions are inconsistent.
Out[1]= {{y[x, t] -> InterpolatingFunction[{{0, 1.}, {0., 1.}}, <>][x, t]}}
In this particular example, a simple way to solve this problem is to replace the initial condition y[x,0]==0 with y[x,0]==If[x>0,0,1] so that it will no longer be inconsistent with the boundary condition. Although this will correct the inconsistency, the use of a discontinuous initial condition will make it difficult to for NDSolve to find a solution. Evidence of this difficulty can be seen in the warning message in the following example, and in the fact that this calculation takes are relatively long time to finish (a few minutes).
In[2]:= NDSolve[{D[y[x, t], t] == D[y[x, t], x, x],
y[x, 0] == If[x > 0, 0, 1], y[0, t] == 1,
Derivative[1, 0][y][1, t] == 0},
y[x, t], {x, 0, 1}, {t, 0, 1}] //Timing
NDSolve::mxsst:
Using maximum number of grid points 200 allowed by the MaxSteps option
Out[2]= {180.371 Second, {{y[x, t] ->
> InterpolatingFunction[{{0, 1.}, {0., 1.}}, <>][x, t]}}}
One way to deal with the problems of this last example is to ignore them. The solution is quite good, and the reason for the warning message is well
understood. You can speed up this calculation by reducing the value of the MaxSteps
option, or by changing the value of the StartingStepSize
option. For example, with StartingStepSize->{.02,Automatic}, the solution will not be quite as good, but the time needed for the calculation will be reduced by about a factor of twelve. The message here is again a reflection of the discontinuity in the initial condition.
In[3]:= NDSolve[{D[y[x, t], t] == D[y[x, t], x, x],
y[x, 0] == If[x > 0, 0, 1], y[0, t] == 1,
Derivative[1, 0][y][1, t] == 0},
y[x, t], {x, 0, 1}, {t, 0, 1},
StartingStepSize -> {.02, Automatic}] //Timing
NDSolve::eerrs:
Warning: Estimated spatial error at t = 0
is much greater than prescribed error tolerance. Grid spacing 0.02
which was based on the value of the StartingStepSize option may be too
large to achieve the desired accuracy or precision.
Out[3]= {14.8878 Second, {{y[x, t] ->
> InterpolatingFunction[{{0, 1.}, {0., 1.}}, <>][x, t]}}}
Another approach to this problem is to use boundary conditions that better approximate the underlying physical problem. In the real world, the quantity represented by y[x,t] probably cannot change discontinuously, either as a function of x or as a function of t, and it may be appropriate to choose boundary conditions which reflect that. The example above might correspond to conduction of heat in a wire, with y[x,t] giving the temperature of the wire at position x and time t. The initial condition y[x,0]==0 indicates that the initial temperature of the wire is 0, and the boundary condition y[0,t]==1 indicates that the temperature at one end of the wire is raised instantaneously from from 0 to 1. Physically, it is not possible for this to happen. In a real wire, non-zero changes in temperature always take place over a finite period of time, and over a finite distance. Inconsistent boundary conditions can also show up as singularities in the solution. Most numerical algorithms, including NDSolve, use smooth functions to approximate the solution, and so will have difficulty in examples where the actual solution is singular. An ideal implementation of this principle would be to incorporate the true physics of the problem into the boundary conditions. In most situations it isn't necessary to be especially careful about this. Almost any smooth function will suffice. Here is an example using y[x,0]==(x-1)^20 in place of y[x,0]==If[x>0,0,1]. The difference between these two expressions can be seen in a plot of the solution, by observing the lack of a discontinuity in the solution for t=0.
In[4]:= sol = NDSolve[{D[y[x, t], t] == D[y[x, t], x, x],
y[x, 0] == (x-1)^20, y[0, t] == 1,
Derivative[1, 0][y][1, t] == 0},
y[x, t], {x, 0, 1}, {t, 0, 1}]
Out[4]= {{y[x, t] -> InterpolatingFunction[{{0, 1.}, {0., 1.}}, <>][x, t]}}
In[5]:= Plot3D[Evaluate[y[x, t] /. sol[[1]]], {x, 0, 1}, {t, 0, 1},
AxesLabel -> {"x", "t", None}]
Out[5]= -SurfaceGraphics-
In general, you may need to experiment with different boundary conditions, and with the options of NDSolve, to get what you want. Despite these difficulties, and the conflict with the underlying physics, there are many applications in which the use of inconsistent boundary and initial conditions is perfectly reasonable, and provides a convenient way to formulate a variety of problems.
Additional Online Documentation:
Mathematica 3.0
http://documents.wolfram.com/v3/RefGuide/NDSolve.html
Mathematica 4.0
http://documents.wolfram.com/v4/RefGuide/NDSolve.html
Questions or comments? Send email to support@wolfram.com.
| |
