Introduction to Advanced Numerical Differential Equation Solving in the Wolfram Language

Overview

The Wolfram Language function (RefLink[NDSolve, paclet:ref/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 (Subscript[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. (RefLink[NDSolve, paclet:ref/NDSolve]) can also solve some differential-algebraic equations (DAEs), which are typically a mix of differential and algebraic equations.

Finding numerical solutions to ordinary differential equations.

(RefLink[NDSolve, paclet:ref/NDSolve]) represents solutions for the functions (Subscript[u, i]) as (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) objects. The (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) objects provide approximations to the (Subscript[u, i]) over the range of values (Subscript[t, min]) to (Subscript[t, max]) for the independent variable t.

In general, (RefLink[NDSolve, paclet:ref/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 (Subscript[t, min]) to (Subscript[t, max]).

In order to get started, (RefLink[NDSolve, paclet:ref/NDSolve]) has to be given appropriate initial or boundary conditions for the (Subscript[u, i]) and their derivatives. These conditions specify values for (Subscript[u, i])[t], and perhaps derivatives (Subscript[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. (RefLink[NDSolve, paclet:ref/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 t==1:

NDSolve[{x'[t]==x[t],x[1]==3},x,{t,0,2}]

When you use (RefLink[NDSolve, paclet:ref/NDSolve]), the initial or boundary conditions you give must be sufficient to determine the solutions for the (Subscript[u, i]) completely. When you use (RefLink[DSolve, paclet:ref/DSolve]) to find symbolic solutions to differential equations, you may specify fewer conditions. The reason is that (RefLink[DSolve, paclet:ref/DSolve]) automatically inserts arbitrary symbolic constants (RefLink[C, paclet:ref/C])[i] to represent degrees of freedom associated with initial conditions that you have not specified explicitly. Since (RefLink[NDSolve, paclet:ref/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 n(^th) derivatives, then you need to either give initial conditions for up to (n-1)(^th) derivatives, or give boundary conditions at n points.

This solves an initial value problem for a second-order equation, which requires two conditions, given at t==0:

NDSolve[{x''[t]==x[t]^2,x[0]==1,x'[0]==0},x,{t,0,2}]

This plots the solution obtained:

Plot[Evaluate[x[t]/.%],{t,0,2}]

Here is a simple boundary value problem:

NDSolve[{y''[x]+x y[x]==0,y[0]==1,y[1]==-1},y,{x,0,1}]

You can use (RefLink[NDSolve, paclet:ref/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:

sol = NDSolve[{x''[t] == y[t] x[t], y'[t] == -(1/(((x[t])^2)+((y[t])^2))), x[0] == 1, x'[0] == 0, y[0] == 0},{x,y}, {t,0,100}]

Here is a plot of both solutions:

Plot[Evaluate[{x[t], y[t]} /. %], {t,0,100}, PlotRange->All, PlotPoints->200]

You can give initial conditions as equations of any kind. If these equations have multiple solutions, (RefLink[NDSolve, paclet:ref/NDSolve]) will generate multiple solutions.

The initial conditions in this case lead to multiple solutions:

NDSolve[{y'[x]^2-y[x]^3==0,y[0]^2==4},y,{x,1}]

(RefLink[NDSolve, paclet:ref/NDSolve]) was not able to find the solution for y'[x]==-(RefLink[Sqrt, paclet:ref/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 (RefLink[NDSolve, paclet:ref/NDSolve]) was able to find. (The upper two solutions are strictly real.)

Plot[Evaluate[Part[Re[y[x] /. %],{1,2,4}]], {x,0, 1}]

(RefLink[NDSolve, paclet:ref/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.

Here is a simple DAE:

NDSolve[{x''[t] + y[t] == x[t], x[t]^2 + y[t]^2 == 1, x[0] == 0, x'[0] == 1},{x,y},{t,0,2}]

Note that while both of the equations have derivative terms, the variable y appears without any derivatives, so (RefLink[NDSolve, paclet:ref/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:

Plot[Evaluate[{x[t], y[t]} /. %],{t,0,1.66}]

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 (Subscript[x, i]).

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

n=25;
(Subscript[x, 0])[t_]:=0;
(Subscript[x, n])[t_]:=1;
eqns = Table[{(Subscript[x, i])'[t] == (n^2)( (Subscript[x, i+1])[t] - 2 (Subscript[x, i])[t] + (Subscript[x, i-1])[t]), (Subscript[x, i])[0] == (((i/n))^10)},{i,n-1}];
vars = Table[(Subscript[x, i])[t],{i,n-1}];
NDSolve[eqns,vars,{t,0,.25}]

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 n.)

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:

ListPlot3D[Table[vars /. First[%], {t,0,.25,.025}]]

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 (RefLink[InterpolatingFunction, paclet:ref/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:

f[x_?VectorQ]:=n^2*ListConvolve[{1,-2,1},x,{2,2},{1,0}];
NDSolve[{X'[t]==f[X[t]],X[0]==(Range[n-1]/n)^10},X,{t,0,.25}]

(RefLink[NDSolve, paclet:ref/NDSolve]) is able to solve some partial differential equations directly when you specify more independent variables.

Finding numerical solutions to partial differential equations.

Here is a solution of the heat equation found directly by (RefLink[NDSolve, paclet:ref/NDSolve]):

NDSolve[{D[u[x,t],t] == D[u[x,t],x,x],u[x,0] == (x^10), u[0,t] == 0, u[1,t] == 1},u,{x,0,1},{t,0,.25}]

Here is a plot of the solution:

Plot3D[Evaluate[First[u[x,t] /. %]],{x,0,1},{t,0,.25}]

(RefLink[NDSolve, paclet:ref/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:

NDSolve[{D[u[t, x, y], t, t] == D[u[t, x, y], x, x] + D[u[t, x, y], y, y] - Sin[u[t, x, y]], u[0, x, y] == Exp[-((x^2)+(y^2))], Derivative[1,0,0][u][0,x,y] == 0, u[t, -5, y] == u[t, 5, y] == 0, u[t, x, -5] == u[t, x, 5] == 0},u,{t,0,3},{x,-5,5},{y,-5,5}]

Here is a plot of the result at t==3:

Plot3D[First[u[3,x,y] /. %],{x,-5,5},{y,-5,5}]

As mentioned earlier, (RefLink[NDSolve, paclet:ref/NDSolve]) works by taking a sequence of steps in the independent variable t. (RefLink[NDSolve, paclet:ref/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 (RefLink[NDSolve, paclet:ref/NDSolve]) will reduce the step size to be able to better track the solution.

(RefLink[NDSolve, paclet:ref/NDSolve]) allows you to specify the precision or accuracy of result you want. In general, (RefLink[NDSolve, paclet:ref/NDSolve]) makes the steps it takes smaller and smaller until the solution reached satisfies either the (RefLink[AccuracyGoal, paclet:ref/AccuracyGoal]) or the (RefLink[PrecisionGoal, paclet:ref/PrecisionGoal]) you give. The setting for (RefLink[AccuracyGoal, paclet:ref/AccuracyGoal]) effectively determines the absolute error to allow in the solution, while the setting for (RefLink[PrecisionGoal, paclet:ref/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 (RefLink[AccuracyGoal, paclet:ref/AccuracyGoal]). By setting (RefLink[AccuracyGoal, paclet:ref/AccuracyGoal])->(RefLink[Infinity, paclet:ref/Infinity]), you tell (RefLink[NDSolve, paclet:ref/NDSolve]) to use (RefLink[PrecisionGoal, paclet:ref/PrecisionGoal]) only. Generally, (RefLink[AccuracyGoal, paclet:ref/AccuracyGoal]) and (RefLink[PrecisionGoal, paclet:ref/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 (RefLink[AccuracyGoal, paclet:ref/AccuracyGoal]) and (RefLink[PrecisionGoal, paclet:ref/PrecisionGoal]).

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

(RefLink[NDSolve, paclet:ref/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 (RefLink[NormFunction, paclet:ref/NormFunction])->f to combine errors in different components. The norm is scaled in terms of the tolerances, given so that (RefLink[NDSolve, paclet:ref/NDSolve]) tries to take steps such that

where (Subscript[err, i]) is the i(^th) component of the error and (Subscript[x, i]) is the i(^th)component of the current solution.

This generates a high-precision solution to a differential equation:

NDSolve[{x'''[t]==x[t],x[0]==1,x'[0]==x''[0]==0},x,{t,1},AccuracyGoal->20,PrecisionGoal->20,WorkingPrecision->25]

Here is the value of the solution at the endpoint:

x[1]/.%

Through its adaptive procedure, (RefLink[NDSolve, paclet:ref/NDSolve]) is able to solve "stiff" differential equations in which there are several components varying with t at extremely different rates.

(RefLink[NDSolve, paclet:ref/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, (RefLink[NDSolve, paclet:ref/NDSolve]) might go on reducing the step size forever, and never terminate. To avoid this problem, the option (RefLink[MaxSteps, paclet:ref/MaxSteps]) specifies the maximum number of steps that (RefLink[NDSolve, paclet:ref/NDSolve]) will ever take in attempting to find a solution. For ordinary differential equations, the default setting is (RefLink[MaxSteps, paclet:ref/MaxSteps])->10000.

(RefLink[NDSolve, paclet:ref/NDSolve]) stops after taking 10,000 steps:

NDSolve[{y'[x]==1/x^2,y[-1]==1},y[x],{x,-1,0}]

There is in fact a singularity in the solution at x=0:

Plot[Evaluate[y[x]/.%],{x,-1,0}]

The default setting (RefLink[MaxSteps, paclet:ref/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 (RefLink[MaxSteps, paclet:ref/MaxSteps]). With the setting (RefLink[MaxSteps, paclet:ref/MaxSteps])->(RefLink[Infinity, paclet:ref/Infinity]) there is no upper limit on the number of steps used.

(RefLink[NDSolve, paclet:ref/NDSolve]) has several different methods built in for computing solutions as well as a mechanism for adding additional methods. With the default setting (RefLink[Method, paclet:ref/Method])->(RefLink[Automatic, paclet:ref/Automatic]), (RefLink[NDSolve, paclet:ref/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 (RefLink[Automatic, paclet:ref/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, (RefLink[NDSolve, paclet:ref/NDSolve]) has a mechanism that allows you to define your own methods and still have the equations and results processed by (RefLink[NDSolve, paclet:ref/NDSolve]) just as for the built-in methods.

When (RefLink[NDSolve, paclet:ref/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 (RefLink[InterpolatingFunction, paclet:ref/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.