Numerical Solution of Differential Equations

The function (RefLink[NDSolve, paclet:ref/NDSolve]) discussed in "Numerical Differential Equations" allows you to find numerical solutions to differential equations. (RefLink[NDSolve, paclet:ref/NDSolve]) handles both single differential equations and sets of simultaneous differential equations. It can handle a wide range of ordinary differential equations as well as some partial differential equations. In a system of ordinary differential equations there can be any number of unknown functions (Subscript[y, i]), but all of these functions must depend on a single "independent variable" x, which is the same for each function. Partial differential equations involve two or more independent variables. (RefLink[NDSolve, paclet:ref/NDSolve]) can also handle differential-algebraic equations that mix differential equations with algebraic ones.

Finding numerical solutions to ordinary differential equations.

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

(RefLink[NDSolve, paclet:ref/NDSolve]) finds solutions iteratively. It starts at a particular value of x, then takes a sequence of steps, trying eventually to cover the whole range (Subscript[x, min]) to (Subscript[x, max]).

In order to get started, (RefLink[NDSolve, paclet:ref/NDSolve]) has to be given appropriate initial or boundary conditions for the (Subscript[y, i]) and their derivatives. These conditions specify values for (Subscript[y, i])[x], and perhaps derivatives (Subscript[y, i])'[x], at particular points x. In general, at least for ordinary differential equations, the conditions you give can be at any x: (RefLink[NDSolve, paclet:ref/NDSolve]) will automatically cover the range (Subscript[x, min]) to (Subscript[x, max]).

This finds a solution for y with x in the range 0 to 2, using an initial condition for y[0]:

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

This still finds a solution with x in the range 0 to 2, but now the initial condition is for y[3]:

NDSolve[{y'[x]==y[x],y[3]==1},y,{x,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}]

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[y, i]) completely. When you use (RefLink[DSolve, paclet:ref/DSolve]) to find symbolic solutions to differential equations, you can get away with specifying fewer initial conditions. The reason is that (RefLink[DSolve, paclet:ref/DSolve]) automatically inserts arbitrary 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 give initial conditions for up to (n-1)(^th) derivatives, or give boundary conditions at n points.

With a third-order equation, you need to give initial conditions for up to second derivatives:

NDSolve[{y'''[x]+8 y''[x]+17 y'[x]+10 y[x]==0,y[0]==6,y'[0]==-20,y''[0]==84},y,{x,0,1}]

This plots the solution obtained:

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

With a third-order equation, you can also give boundary conditions at three points:

NDSolve[{y'''[x]+Sin[x]==0,y[0]==4,y[1]==7,y[2]==0},y,{x,0,2}]

The Wolfram Language allows you to use any appropriate linear combination of function values and derivatives as boundary conditions:

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

In most cases, all the initial conditions you give must involve the same value of x, say (Subscript[x, 0]). As a result, you can avoid giving both (Subscript[x, min]) and (Subscript[x, max]) explicitly. If you specify your range of x as {x,(Subscript[x, 1])}, then the Wolfram Language will automatically generate a solution over the range (Subscript[x, 0]) to (Subscript[x, 1]).

This generates a solution over the range 0 to 2:

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

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

The initial conditions in this case lead to multiple solutions:

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

Here is a plot of all the solutions:

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

You can use (RefLink[NDSolve, paclet:ref/NDSolve]) to solve systems of coupled differential equations.

This finds a numerical solution to a pair of coupled equations:

sol=NDSolve[{x'[t]==-y[t]-x[t]^2,y'[t]==2 x[t]-y[t],x[0]==y[0]==1},{x,y},{t,10}]

This plots the solution for y from these equations:

Plot[Evaluate[y[t]/.sol],{t,0,10},PlotRange->All]

This generates a parametric plot using both x and y:

ParametricPlot[Evaluate[{x[t],y[t]}/.sol],{t,0,10},PlotRange->All]

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

This constructs a set of five coupled differential equations and initial conditions:

eqns=Join[Table[y[i]'[x]==y[i-1][x]-y[i][x],{i,2,4}],{y[1]'[x]==-y[1][x],y[5]'[x]==y[4][x],y[1][0]==1},Table[y[i][0]==0,{i,2,5}]]

This solves the equations:

NDSolve[eqns,Table[y[i],{i,5}],{x,10}]

Here is a plot of the solutions:

Plot[Evaluate[Table[y[i][x],{i,5}]/.%],{x,0,10}]

(RefLink[NDSolve, paclet:ref/NDSolve]) can handle functions whose values are lists or arrays. If you give initial conditions like y[0]=={(Subscript[v, 1]),(Subscript[v, 2]),…,(Subscript[v, n])}, then (RefLink[NDSolve, paclet:ref/NDSolve]) will assume that y is a function whose values are lists of length n.

This solves a system of four coupled differential equations:

NDSolve[{y''[x]==-RandomReal[{0,1},{4,4}].y[x],y[0]==y'[0]==Table[1,{4}]},y,{x,0,8}]

Here are the solutions:

With[{s=y[x]/.First[%]},Plot[{s[[1]],s[[2]],s[[3]],s[[4]]},{x,0,8},PlotRange->All]]

Special options for (RefLink[NDSolve, paclet:ref/NDSolve]).

(RefLink[NDSolve, paclet:ref/NDSolve]) has many methods for solving equations, but essentially all of them at some level work by taking a sequence of steps in the independent variable x, and using 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 or change the method so as to be able to track the solution better.

This solves a differential equation in which the derivative has a discontinuity:

NDSolve[{y'[x]==If[x<0,1/(x-1),1/(x+1)],y[-5]==5},y,{x,-5,5}]

(RefLink[NDSolve, paclet:ref/NDSolve]) reduced the step size around x==0 so as to reproduce the kink accurately:

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

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

In these equations, y varies much more rapidly than z:

sol=NDSolve[{y'[x]==-40 y[x],z'[x]==-z[x]/10,y[0]==z[0]==1},{y,z},{x,0,1}]

(RefLink[NDSolve, paclet:ref/NDSolve]) nevertheless tracks both components successfully:

Plot[Evaluate[{y[x],z[x]}/.sol],{x,0,1},PlotRange->All]

(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 or the integration interval is too big. For the first case, (RefLink[NDSolve, paclet:ref/NDSolve]) limits the smallest step size that it will consider as significant for a given integration interval. For the second case, 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])->(RefLink[Automatic, paclet:ref/Automatic]). With the (RefLink[Automatic, paclet:ref/Automatic]) setting, (RefLink[NDSolve, paclet:ref/NDSolve]) will estimate how many steps are needed to solve the equation at hand based on the initial step sizes taken.

(RefLink[NDSolve, paclet:ref/NDSolve]) stops when the step size becomes too small:

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 for (RefLink[MaxSteps, paclet:ref/MaxSteps]) should be sufficient for most equations with smooth solutions. When solutions have a complicated structure, however, you may occasionally 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.

To take the solution to the Lorenz equations this far, you need to remove the default bound on (RefLink[MaxSteps, paclet:ref/MaxSteps]):

NDSolve[{x'[t]==-3 (x[t]-y[t]),y'[t]==-x[t] z[t]+26.5 x[t]-y[t],z'[t]==x[t] y[t]-z[t],x[0]==z[0]==0,y[0]==1},{x,y,z},{t,0,200},MaxSteps->Infinity]

Here is a parametric plot of the solution in three dimensions:

ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.%],{t,0,200},PlotPoints->10000,ColorFunction->(ColorData["Rainbow"][#4]&)]

When (RefLink[NDSolve, paclet:ref/NDSolve]) solves a particular set of differential equations, it always tries to choose a step size appropriate for those equations. In some cases, the very first step that (RefLink[NDSolve, paclet:ref/NDSolve]) makes may be too large, and it may miss an important feature in the solution. To avoid this problem, you can explicitly set the option (RefLink[StartingStepSize, paclet:ref/StartingStepSize]) to specify the size to use for the first step.

The equations you give to (RefLink[NDSolve, paclet:ref/NDSolve]) do not necessarily all have to involve derivatives; they can also just be algebraic. You can use (RefLink[NDSolve, paclet:ref/NDSolve]) to solve many such differential-algebraic equations.

This solves a system of differential-algebraic equations:

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

Here is the solution:

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

Finding numerical solutions to partial differential equations.

This finds a numerical solution to the wave equation. The result is a two-dimensional interpolating function:

NDSolve[{D[u[t,x],t,t]==D[u[t,x],x,x],u[0,x]==Exp[-x^2],Derivative[1,0][u][0,x]==0,u[t,-6]==u[t,6]},u,{t,0,6},{x,-6,6}]

This generates a plot of the result:

Plot3D[Evaluate[u[t,x]/.First[%]],{t,0,6},{x,-6,6},PlotPoints->50]

This finds a numerical solution to a nonlinear wave equation:

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

Here is a 3D plot of the result:

Plot3D[Evaluate[u[t,x]/.First[%]],{t,0,10},{x,-10,10}]

This is a higher-resolution density plot of the solution:

DensityPlot[Evaluate[u[10-t,x]/.First[%%]],{x,-10,10},{t,0,10},PlotPoints->200,Mesh->False]

Here is a version of the equation in 2+1 dimensions:

eqn=D[u[t,x,y],t,t]==D[u[t,x,y],x,x]+D[u[t,x,y],y,y]/2+(1-u[t,x,y]^2)(1+2u[t,x,y])

This solves the equation:

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

This generates a list of plots of the solution:

Table[Plot3D[Evaluate[u[t,x,y]/.First[%]],{x,-5,5},{y,-5,5},PlotRange->All,PlotPoints->100,Mesh->False],{t,1,4}]