This is documentation for Mathematica 3, which was
based on an earlier version of the Wolfram Language.
View current documentation (Version 11.1)
 Documentation / Mathematica / Built-in Functions / New in Version 3.0 / Numerical Computation  /

  • 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 , , , ... , x , xmin , xmax ] finds numerical solutions for the functions¬†.
  • NDSolve gives results in terms of InterpolatingFunction objects.
  • NDSolve[ eqns , y [ x ], x , xmin , xmax ] gives solutions for y [ x ] rather than for the function y itself.
  • Differential equations must be stated in terms of derivatives such as y '[ x ], obtained with D, not total derivatives obtained with Dt.
  • NDSolve solves a wide range of ordinary differential equations, and some partial differential equations.
  • In ordinary differential equations the functions must depend only on the single variable x. In partial differential equations they may depend on more than one variable.
  • The differential equations must contain enough initial or boundary conditions to determine the solutions for the completely.
  • Initial and boundary conditions are typically stated in form y [ ] == , y '[ ] == , etc., but may consist of more complicated equations.
  • The point that appears in the initial or boundary conditions need not lie in the range xmin to xmax over which the solution is sought.
  • The differential equations in NDSolve can involve complex numbers.
  • The following options can be given:
  • The default setting for MaxSteps is 1000 for ordinary differential equations, and 200 for partial differential equations.
  • NDSolve stops when either the AccuracyGoal or the PrecisionGoal specified is met.
  • The default setting of Automatic for AccuracyGoal and PrecisionGoal yields goals equal to the setting for WorkingPrecision minus 10 digits.
  • AccuracyGoal effectively specifies the absolute error allowed in solutions, while PrecisionGoal specifies the relative error.
  • If solutions must be followed accurately when their values are close to zero, AccuracyGoal should be set larger, or to Infinity.
  • The setting for InterpolationPrecision specifies the number of digits of precision to use inside the InterpolatingFunction object generated by NDSolve.
  • The default setting of Automatic for InterpolationPrecision uses the current setting for WorkingPrecision.
  • See the Mathematica book: Section 3.9.1,¬†Section 3.9.7.
  • See also Implementation NotesA.9.44.21MainBookLinkOldButtonDataA.9.44.21.
  • See also: DSolve, NIntegrate.

    Further Examples

    Ordinary differential equations: Basic usage
    This command finds a numerical approximation to a function that is equal to its first derivative at each point x between and , and that has the value when x is .
    NDSolve returns a rule to replace y by an InterpolatingFunction object.



    Mathematica expresses the closed-form solution of this interesting differential equation as a pure function.



    We rewrite it as an ordinary algebraic expression. The solution has a singularity at .



    Now we solve the same equation approximately, using NDSolve instead of DSolve. As we can see from the closed-form solution, the warning is appropriate.


    NDSolve::ndsz: At x == 0.999982, step size is effectively zero; singularity suspected.


    We can use an InterpolatingFunction object like any other function that evaluates to a number. Here are three ways to check the boundary condition.







    The solution was not found all the way to , but the InterpolatingFunction object still allows you to evaluate with arguments outside its range. You can see that the approximate value is quite large, though smaller than the correct value of infinity!


    InterpolatingFunction::dmval: Input value {1} lies outside the range of data in the interpolating function. Extrapolation will be used.


    It is also easy to make plots of solutions. Here is how to see a graph of an approximation to the function which is the reciprocal of its derivative. The Evaluate command saves time by substituting the InterpolatingFunction once, instead of for each number used to generate the plot.
    Not too surprisingly, this plot looks very much like the square root function. What do you think would happen if you tried to solve in the other direction?

    Evaluate the cell to see the graphic.


    Mathematica can handle higher order equations. Here is a plot of the interesting solution of a third order equation.
    Note that you have to specify enough initial conditions to uniquely determine the solution.

    Evaluate the cell to see the graphic.


    Solving systems of equations works similarly. For systems of two equations a so-called phase plot is often a good way to visualize the solution. Here is a phase plot that describes the motion of a weakly damped pendulum.

    Evaluate the cell to see the graphic.


    Ordinary differential equations: Options
    This way of solving the differential equation for the cosine function doesn't work over such a long interval.


    NDSolve::mxst: Maximum number of 1000 steps reached at the point x == 199.224.


    You need to increase the setting of the MaxSteps option. As expected, at the end of the interval the solution is close to .



    The Method option allows you to specify which method NDSolve uses to approximate the solution. The default method, Automatic, automatically switches between Gear and Adams, depending on stiffness. Another possibility for equations which are not stiff is RungeKutta. For some problems, the Runge-Kutta method can find the solution using fewer steps.
    For the Rössler equations, the Runge-Kutta method needs about half as many steps as the default method.

    Evaluate the cell to see the graphic.


    One way to get a very precise solution of an ODE is to give a sufficiently high value for the WorkingPrecision option. Note that AccuracyGoal and PrecisionGoal default to 10 less than the value of WorkingPrecision when it is greater than $MachinePrecision.
    You need to be careful not to use too high a value for WorkingPrecision, because as working precision increases, not only does the time taken for each arithmetic operation and function evaluation increase, but the number of steps typically increases exponentially also.
    These commands compare a known exact solution with solutions computed with different values of WorkingPrecision.




    The values of x at which the steps are taken is kept in the third part of the InterpolatingFunction object. This is why the Length command above gives the number of steps.



    Ordinary differential equations: Boundary value problems
    Simple linear boundary value problems can be solved by constructing the boundary value equations appropriately. If the order of the equation is you need to know the values of some combination of the function or its derivatives at points.
    This normalized equation describes the effect of a wave incident on the right edge of an optical medium (where x is ).

    Evaluate the cell to see the graphic.


    Here is a third-order equation where the function values and combinations of the derivatives are known at three points.

    Evaluate the cell to see the graphic.


    Not all linear equations with boundary values can be solved by the method that is implemented. Nonlinear equations cannot be solved, either.


    NDSolve::unsol: Not possible to initiate boundary value problem with the chasing method

    NDSolve::unsol: Not possible to initiate boundary value problem with the chasing method

    Not possible to initiate boundary value problem
    with the chasing method

    Not possible to initiate boundary value problem
    with the chasing method

    Further output of NDSolve::unsol
    will be suppressed during this calculation.


    Partial differential equations: Basic usage
    For the first three examples, consider the one-dimensional heat equation. With fixed boundary conditions this is a model for diffusion of heat in an insulated rod with the temperatures at the endpoints held fixed.
    This command solves the heat equation with the left end () held at fixed temperature , the right end () held at fixed temperature , and an initial heat profile given by the quadratic in x.



    Just as for ODE's, the solution is returned as a rule to replace the dependent variable with an InterpolatingFunction object. The only noticeable difference is that the InterpolatingFunction object is two-dimensional. That is to say, it takes two arguments to evaluate, with the arguments in the same order as the variables in the NDSolve command.
    The boundary conditions can be a linear combination of Dirichlet and Neumann type conditions. For example, this command solves a model for the diffusion of heat in a one-dimensional insulated rod with the left end held at a constant temperature and the right end radiating into free space. The Neumann boundary condition was entered using the Derivative operator. An equivalent way to give the condition is ((u[x,t] + D[u[x,t],x]) /. x->1) == 0.

    Evaluate the cell to see the graphic.


    An even more interesting case is when the boundary conditions are time-dependent. For example, entering these command produces a plot of a solution with the temperature at the left edge varying sinusoidally.

    Evaluate the cell to see the graphic.


    It is also possible to compute the solution of some systems of PDE's. These commands solve a system of mixed parabolic-hyperbolic type and produce separate contour plots for each of the dependent variables.

    Evaluate the cell to see the graphics.


    Periodic boundary conditions are frequently convenient for numerical solutions. You can tell NDSolve to solve with periodic boundary conditions by specifying that the values of the solution are to be equal at the left and right edges of the domain in one independent variable for all values of the other independent variable.
    For example, an interesting solution with periodic solutions is the nonlinear Schrödinger equation. These commands set up a periodic initial condition (it happens to be a soliton), and computes the solution, and then produces plots showing the modulus and real and imaginary parts of the solution.

    Evaluate the cell to see the graphics.


    Partial differential equations: Limitations
    NDSolve uses a method for solving PDE's that is called the numerical method of lines. It discretizes in one variable to make a system of ODE's. This system is then solved using the ODE methods built into NDSolve.
    For the method to work, an initial function must be specified for one variable and boundary values may be specified for the other variable. The initial function is used to find the initial conditions for the system of the ODE's. Boundary and initial values may be specified on at most three sides of a rectangle.
    The method has the advantage that it can solve a reasonably large class of equations. However, there are types of equations which it cannot solve.
    Because elliptic problems are ill-posed unless boundary values are specified on all sides of a region, this method cannot find solutions of elliptic problems. A classic example is Laplace's equation. Entering this command indicates what will happen when you try to do something of this sort. Typically ill-posedness will appear as numerical instability, and looking at the scale on the plot produced indicates that you should heed the message seriously!

    Evaluate the cell to see the graphic.


    Another class of problems that the method cannot currently handle are those that form singularities in the solution. The discretization is done with finite difference methods, so fronts may be incorrectly resolved or completely lost. For example, Bergers' equation is a model for some of the features of gas dynamics, including shock formation. Entering these commands will produce a plot that show the oscillations typical of the interaction between finite difference methods and fronts.

    Evaluate the cell to see the graphic.


    Partial differential equations: Options
    Many of the options that control the ODE solver also apply to the PDE solver, though often in a different way. The solutions for PDE's are computed in two stages. First the equation is discretized, and then the resulting system of ODE's is solved using NDSolve's built-in method. If you want different option values for the two stages, you can specify the option value as a list. (The order of the independent variables in the command determines to what variable the options apply.)
    The following options to NDSolve can be used in such a list: AccuracyGoal, PrecisionGoal, DifferenceOrder, MaxSteps, MaxStepSize, StartingStepSize.
    In the discretization stage, the default used is fourth-order finite differences. In some cases fourth-order is not optimal; you can control this with the DifferenceOrder option. When the equation has high spatial differential order (eg. Airy's equation), it is better to increase the order.
    Entering these commands computes a solution to Airy's equation with periodic boundary conditions and produces plots that you can view as an animation. Since the independent variable is the one with an initial function, it is the one for which we specify the difference order of . (The variable x appears after t in the list of arguments of NDSolve.)

    Evaluate the cell to see the graphic.


    We clean up by clearing the definitions that were made.