Introduction to Advanced Numerical Differential Equation Solving in the Wolfram Language

Overview

The Wolfram Language function NDSolve is a general numerical differential equation solver. It can handle a wide range of ordinary differential equations (ODEs) as well as some partial differential equations (PDEs). In a system of ordinary differential equations there can be any number of unknown functions u i , but all of these functions must depend on a single "independent variable" t , which is the same for each function. Partial differential equations involve two or more independent variables. NDSolve can also solve some differential-algebraic equations (DAEs), which are typically a mix of differential and algebraic equations.

NDSolve[ { eqn 1 , eqn 2 , … } , u , { t , t min , t max } ] find a numerical solution for the function u with t in the range t min to t max
NDSolve[ { eqn 1 , eqn 2 , … } , { u 1 , u 2 , … } , { t , t min , t max } ] find numerical solutions for several functions u i

Finding numerical solutions to ordinary differential equations.

NDSolve represents solutions for the functions u i as InterpolatingFunction objects. The InterpolatingFunction objects provide approximations to the u i over the range of values t min to t max for the independent variable t .

In general, NDSolve finds solutions iteratively. It starts at a particular value of t , then takes a sequence of steps, trying eventually to cover the whole range t min to t max .

In order to get started, NDSolve has to be given appropriate initial or boundary conditions for the u i and their derivatives. These conditions specify values for u i [ t ] , and perhaps derivatives u i '[ t ] , at particular points t . When there is only one t at which conditions are given, the equations and initial conditions are collectively referred to as an initial value problem. A boundary value occurs when there are multiple points t . NDSolve can solve nearly all initial value problems that can symbolically be put in normal form (i.e. are solvable for the highest derivative order), but only linear boundary value problems.

This finds a solution for x with t in the range 0 to 2, using an initial condition for x at t1:

When you use NDSolve, the initial or boundary conditions you give must be sufficient to determine the solutions for the u i completely. When you use DSolve to find symbolic solutions to differential equations, you may specify fewer conditions. The reason is that DSolve automatically inserts arbitrary symbolic constants C[ i ] to represent degrees of freedom associated with initial conditions that you have not specified explicitly. Since NDSolve must give a numerical solution, it cannot represent these kinds of additional degrees of freedom. As a result, you must explicitly give all the initial or boundary conditions that are needed to determine the solution.

In a typical case, if you have differential equations with up to th derivatives, then you need to either give initial conditions for up to th derivatives, or give boundary conditions at points.

This solves an initial value problem for a second-order equation, which requires two conditions, given at t0:
This plots the solution obtained:
Here is a simple boundary value problem:

You can use NDSolve to solve systems of coupled differential equations as long as each variable has the appropriate number of conditions.

This finds a numerical solution to a pair of coupled equations:
Here is a plot of both solutions:

You can give initial conditions as equations of any kind. If these equations have multiple solutions, NDSolve will generate multiple solutions.

The initial conditions in this case lead to multiple solutions:

NDSolve was not able to find the solution for y'[x] -Sqrt[y[x]^3] , y[0] -2 because of problems with the branch cut in the square root function.

This shows the real part of the solutions that NDSolve was able to find. (The upper two solutions are strictly real.)

NDSolve can solve a mixed system of differential and algebraic equations, referred to as differential-algebraic equations (DAEs). In fact, the example given is a sort of DAE, where the equations are not expressed explicitly in terms of the derivatives. Typically, however, in DAEs, you are not able to solve for the derivatives at all and the problem must be solved using a different method entirely.

Note that while both of the equations have derivative terms, the variable y appears without any derivatives, so NDSolve issues a warning message. When the usual substitution to convert to first-order equations is made, one of the equations does indeed become effectively algebraic.

Also, since y only appears algebraically, it is not necessary to give an initial condition to determine its values. Finding initial conditions that are consistent with DAEs can, in fact, be quite difficult. The tutorial "Numerical Solution of Differential-Algebraic Equations" has more information.

This shows a plot of the solutions:

From the plot, you can see that the derivative of y is tending to vary arbitrarily fast. Even though it does not explicitly appear in the equations, this condition means that the solver cannot continue further.

Unknown functions in differential equations do not necessarily have to be represented by single symbols. If you have a large number of unknown functions, for example, you will often find it more convenient to give the functions names like x [ i ] or x i .

This constructs a set of 25 coupled differential equations and initial conditions and solves them:

This actually computes an approximate solution of the heat equation for a rod with constant temperatures at either end of the rod. (For more accurate solutions, you can increase .)

The result is an approximate solution to the heat equation for a 1-dimensional rod of length 1 with constant temperature maintained at either end. This shows the solutions considered as spatial values as a function of time:

An unknown function can also be specified to have a vector (or matrix) value. The dimensionality of an unknown function is taken from its initial condition. You can mix scalar and vector unknown functions as long as the equations have consistent dimensionality according to the rules of Wolfram Language arithmetic. The InterpolatingFunction result will give values with the same dimensionality as the unknown function. Using nonscalar variables is very convenient when a system of differential equations is governed by a process that may be difficult or inefficient to express symbolically.

This uses a vector-valued unknown function to solve the same system as earlier:

NDSolve is able to solve some partial differential equations directly when you specify more independent variables.

NDSolve[ { eqn 1 , eqn 2 , … } , u , { t , t min , t max } , { x , x min , x max } , … ]
solve a system of partial differential equations for a function u [ t , x … ] with t in the range t min to t max and x in the range x min to x max , …
NDSolve[ { eqn 1 , eqn 2 , … } , { u 1 , u 2 , … } , { t , t min , t max } , { x , x min , x max } , … ]
solve a system of partial differential equations for several functions u i

Finding numerical solutions to partial differential equations.

Here is a solution of the heat equation found directly by NDSolve:
Here is a plot of the solution:

NDSolve currently uses the numerical method of lines to compute solutions to partial differential equations. The method is restricted to problems that can be posed with an initial condition in at least one independent variable. For example, the method cannot solve elliptic PDEs such as Laplace's equation because these require boundary values. For the problems it does solve, the method of lines is quite general, handling systems of PDEs or nonlinearity well, and often quite fast. Details of the method are given in "The Numerical Method of Lines".

This finds a numerical solution to a generalization of the nonlinear sine‐Gordon equation to two spatial dimensions with periodic boundary conditions:
Here is a plot of the result at t3:

As mentioned earlier, NDSolve works by taking a sequence of steps in the independent variable t . NDSolve uses an adaptive procedure to determine the size of these steps. In general, if the solution appears to be varying rapidly in a particular region, then NDSolve will reduce the step size to be able to better track the solution.

NDSolve allows you to specify the precision or accuracy of result you want. In general, NDSolve makes the steps it takes smaller and smaller until the solution reached satisfies either the AccuracyGoal or the PrecisionGoal you give. The setting for AccuracyGoal effectively determines the absolute error to allow in the solution, while the setting for PrecisionGoal determines the relative error. If you need to track a solution whose value comes close to zero, then you will typically need to increase the setting for AccuracyGoal. By setting AccuracyGoal->Infinity, you tell NDSolve to use PrecisionGoal only. Generally, AccuracyGoal and PrecisionGoal are used to control the error local to a particular time step. For some differential equations, this error can accumulate, so it is possible that the precision or accuracy of the result at the end of the time interval may be much less than what you might expect from the settings of AccuracyGoal and PrecisionGoal.

NDSolve uses the setting you give for WorkingPrecision to determine the precision to use in its internal computations. If you specify large values for AccuracyGoal or PrecisionGoal, then you typically need to give a somewhat larger value for WorkingPrecision. With the default setting of Automatic, both AccuracyGoal and PrecisionGoal are equal to half of the setting for WorkingPrecision.

NDSolve uses error estimates for determining whether it is meeting the specified tolerances. When working with systems of equations, it uses the setting of the option NormFunction->f to combine errors in different components. The norm is scaled in terms of the tolerances, given so that NDSolve tries to take steps such that

where is the th component of the error and is the thcomponent of the current solution.

This generates a high‐precision solution to a differential equation:
Here is the value of the solution at the endpoint:

Through its adaptive procedure, NDSolve is able to solve "stiff" differential equations in which there are several components varying with t at extremely different rates.

NDSolve follows the general procedure of reducing step size until it tracks solutions accurately. There is a problem, however, when the true solution has a singularity. In this case, NDSolve might go on reducing the step size forever, and never terminate. To avoid this problem, the option MaxSteps specifies the maximum number of steps that NDSolve will ever take in attempting to find a solution. For ordinary differential equations, the default setting is MaxSteps->10000.

NDSolve stops after taking 10,000 steps:
There is in fact a singularity in the solution at x =0:

The default setting MaxSteps->10000 should be sufficient for most equations with smooth solutions. When solutions have a complicated structure, however, you may sometimes have to choose larger settings for MaxSteps. With the setting MaxSteps->Infinity there is no upper limit on the number of steps used.

NDSolve has several different methods built in for computing solutions as well as a mechanism for adding additional methods. With the default setting Method->Automatic, NDSolve will choose a method which should be appropriate for the differential equations. For example, if the equations have stiffness, implicit methods will be used as needed, or if the equations make a DAE, a special DAE method will be used. In general, it is not possible to determine the nature of solutions to differential equations without actually solving them: thus, the default Automatic methods are good for solving as wide variety of problems, but the one chosen may not be the best one available for your particular problem. Also, you may want to choose methods, such as symplectic integrators, which preserve certain properties of the solution.

Choosing an appropriate method for a particular system can be quite difficult. To complicate it further, many methods have their own settings, which can greatly affect solution efficiency and accuracy. Much of this documentation consists of descriptions of methods to give you an idea of when they should be used and how to adjust them to solve particular problems. Furthermore, NDSolve has a mechanism that allows you to define your own methods and still have the equations and results processed by NDSolve just as for the built-in methods.

When NDSolve computes a solution, there are typically three phases. First, the equations are processed, usually into a function that represents the right-hand side of the equations in normal form. Next, the function is used to iterate the solution from the initial conditions. Finally, data saved during the iteration procedure is processed into one or more InterpolatingFunction objects. Using functions in the NDSolve` context, you can run these steps separately and, more importantly, have more control over the iteration process. The steps are tied by an NDSolve`StateData object, which keeps all of the data necessary for solving the differential equations.