({{, (RefLink[NDSolve, paclet:ref/NDSolve])[eqns,u,{x,(Subscript[x, min]),(Subscript[x, max])}]
finds a numerical solution to the ordinary differential equations eqns for the function u with the independent variable x in the range (Subscript[x, min]) to (Subscript[x, max]). }, {, (RefLink[NDSolve, paclet:ref/NDSolve])[eqns,u,{x,(Subscript[x, min]),(Subscript[x, max])},{y,(Subscript[y, min]),(Subscript[y, max])}]
solves the partial differential equations eqns over a rectangular region.}, {, (RefLink[NDSolve, paclet:ref/NDSolve])[eqns,u,{x,y}∈Ω]
solves the partial differential equations eqns over the region Ω.}, {, (RefLink[NDSolve, paclet:ref/NDSolve])[eqns,u,{t,(Subscript[t, min]),(Subscript[t, max])},{x,y}∈Ω]
solves the time-dependent partial differential equations eqns over the region Ω.}, {, (RefLink[NDSolve, paclet:ref/NDSolve])[eqns,{(Subscript[u, 1]),(Subscript[u, 2]),…},…]
solves for the functions (Subscript[u, i]). }})

OpenersNotesSection(Spacer1)Details and Options

(RefLink[NDSolve, paclet:ref/NDSolve]) gives results in terms of (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) objects.

(RefLink[NDSolve, paclet:ref/NDSolve])[eqns,u[x],{x,(Subscript[x, min]),(Subscript[x, max])}] gives solutions for u[x] rather than for the function u itself.

Differential equations must be stated in terms of derivatives such as u'[x], obtained with (RefLink[D, paclet:ref/D]), not total derivatives obtained with (RefLink[Dt, paclet:ref/Dt]).

Partial differential equations may also be specified using the differential operators (RefLink[Grad, paclet:ref/Grad]) (∇), (RefLink[Div, paclet:ref/Div]) (∇.), (RefLink[Laplacian, paclet:ref/Laplacian]) ((∇^2)), and (RefLink[Curl, paclet:ref/Curl]) (∇×). Typically these operators are used as in (RefLink[Inactive, paclet:ref/Inactive])[op] to keep the operator form from evaluating.

(RefLink[NDSolve, paclet:ref/NDSolve]) solves a wide range of ordinary differential equations as well as many partial differential equations.

(RefLink[NDSolve, paclet:ref/NDSolve]) can also solve many delay differential equations.

In ordinary differential equations, the functions (Subscript[u, i]) must depend only on the single variable t. In partial differential equations, they may depend on more than one variable.

(RefLink[WhenEvent, paclet:ref/WhenEvent])[event,action] may be included in the equations eqns to specify an action that occurs when event becomes (RefLink[True, paclet:ref/True]).

The differential equations must contain enough initial or boundary conditions to determine the solutions for the (Subscript[u, i]) completely.

Initial and boundary conditions are typically stated in the form u[(Subscript[x, 0])]==(Subscript[c, 0]), u'[(Subscript[x, 0])]==(Subscript[dc, 0]), etc., but may consist of more complicated equations.

The (Subscript[c, 0]), (Subscript[dc, 0]), etc. can be lists, specifying that u[x] is a function with vector or general list values.

Periodic boundary conditions can be specified using u[(Subscript[x, 0])]==u[(Subscript[x, 1])].

The point (Subscript[x, 0]) that appears in the initial or boundary conditions need not lie in the range (Subscript[x, min]) to (Subscript[x, max]) over which the solution is sought.

Boundary values may also be specified using (RefLink[DirichletCondition, paclet:ref/DirichletCondition]) and (RefLink[NeumannValue, paclet:ref/NeumannValue]).

In delay differential equations, initial history functions are given in the form u[x/;x<(Subscript[x, 0])]==(Subscript[c, 0]), where (Subscript[c, 0]) is in general a function of x.

The differential equations in (RefLink[NDSolve, paclet:ref/NDSolve]) can involve complex numbers.

(RefLink[NDSolve, paclet:ref/NDSolve]) can solve many differential-algebraic equations, in which some of the eqns are purely algebraic, or some of the variables are implicitly algebraic.

The (Subscript[u, i]) can be functions of the dependent variables and need not include all such variables.

The following options can be given:

(RefLink[NDSolve, paclet:ref/NDSolve]) adapts its step size so that the estimated error in the solution is just within the tolerances specified by (RefLink[PrecisionGoal, paclet:ref/PrecisionGoal]) and (RefLink[AccuracyGoal, paclet:ref/AccuracyGoal]).

The option (RefLink[NormFunction, paclet:ref/NormFunction])->f specifies that the estimated errors for each of the (Subscript[u, i]) should be combined using f[{(Subscript[e, 1]),(Subscript[e, 2]),…}].

(RefLink[AccuracyGoal, paclet:ref/AccuracyGoal]) effectively specifies the absolute local error allowed at each step in finding a solution, while (RefLink[PrecisionGoal, paclet:ref/PrecisionGoal]) specifies the relative local error.

If solutions must be followed accurately when their values are close to zero, (RefLink[AccuracyGoal, paclet:ref/AccuracyGoal]) should be set larger, or to (RefLink[Infinity, paclet:ref/Infinity]).

The default setting of (RefLink[Automatic, paclet:ref/Automatic]) for (RefLink[AccuracyGoal, paclet:ref/AccuracyGoal]) and (RefLink[PrecisionGoal, paclet:ref/PrecisionGoal]) is equivalent to (RefLink[WorkingPrecision, paclet:ref/WorkingPrecision])/2.

The default setting of (RefLink[Automatic, paclet:ref/Automatic]) for (RefLink[MaxSteps, paclet:ref/MaxSteps]) estimates the maximum number of steps to be taken by (RefLink[NDSolve, paclet:ref/NDSolve]), depending on start and stop time and an estimate of the step size. Should this not be possible, a fixed number of steps is taken.

The setting for (RefLink[MaxStepFraction, paclet:ref/MaxStepFraction]) specifies the maximum step to be taken by (RefLink[NDSolve, paclet:ref/NDSolve]) as a fraction of the range of values for each independent variable.

With (RefLink[DependentVariables, paclet:ref/DependentVariables])->(RefLink[Automatic, paclet:ref/Automatic]), (RefLink[NDSolve, paclet:ref/NDSolve]) attempts to determine the dependent variables by analyzing the equations given.

(RefLink[NDSolve, paclet:ref/NDSolve]) typically solves differential equations by going through several different stages, depending on the type of equations. With (RefLink[Method, paclet:ref/Method])->{(Subscript[s, 1])->(Subscript[m, 1]),(Subscript[s, 2])->(Subscript[m, 2]),…}, stage  (Subscript[s, i]) is handled by method (Subscript[m, i]). The actual stages used and their order are determined by (RefLink[NDSolve, paclet:ref/NDSolve]), based on the problem to solve.

Possible solution stages are:

With (RefLink[Method, paclet:ref/Method])->(Subscript[m, 1]) or (RefLink[Method, paclet:ref/Method])->{(Subscript[m, 1]),(Subscript[s, 2])->(Subscript[m, 2]),…}, the method (Subscript[m, 1]) is assumed to be for time integration, so (RefLink[Method, paclet:ref/Method])->(Subscript[m, 1]) is equivalent to (RefLink[Method, paclet:ref/Method])->{"TimeIntegration"->(Subscript[m, 1])}.

Possible explicit time integration settings for the (RefLink[Method, paclet:ref/Method]) option include:

With (RefLink[Method, paclet:ref/Method])->{"controller",(RefLink[Method, paclet:ref/Method])->"submethod"} or (RefLink[Method, paclet:ref/Method])->{"controller",(RefLink[Method, paclet:ref/Method])->{(Subscript[m, 1]),(Subscript[m, 2]),…}}, possible controller methods include:

Methods used mainly as submethods include:

The setting (RefLink[InterpolationOrder, paclet:ref/InterpolationOrder])->(RefLink[All, paclet:ref/All]) specifies that (RefLink[NDSolve, paclet:ref/NDSolve]) should generate solutions that use interpolation of the same order as the underlying method used. »

OpenersPrimaryExamplesSection(Spacer1)Examples  (111)

OpenersExampleSection0(Spacer1)Basic Examples  (7)

Solve a first-order ordinary differential equation:

s=NDSolve[{y'[x]==y[x]Cos[x+y[x]],y[0]==1},y,{x,0,30}]

Use the solution in a plot:

Plot[Evaluate[y[x]/.s],{x,0,30},PlotRange->All]

Use the function and its derivative in a plot:

ParametricPlot[Evaluate[{y[x],y'[x]}/.s],{x,0,20}]

Find specific values:

{y[10.5],y'[12.5]}/.s

Second-order nonlinear ordinary differential equation:

s=NDSolve[{y''[x]+ Sin[y[x]]y[x]==0,y[0]==1,y'[0]==0},y,{x,0,30}]

Plot the function and its first two derivatives:

Plot[Evaluate[{y[x],y'[x],y''[x]}/.s],{x,0,30},PlotStyle->Automatic]

System of ordinary differential equations:

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

ParametricPlot[Evaluate[{x[t],y[t]}/.s],{t,0,20}]

Solve the heat equation in one dimension:

NDSolve[{D[u[t,x],t]==D[u[t,x],x,x],u[0,x]==0,u[t,0]==Sin[t],u[t,5]==0},u,{t,0,10},{x,0,5}]

Plot3D[Evaluate[u[t,x]/.%],{t,0,10},{x,0,5}, PlotRange->All]

Alternative form of equation:

NDSolve[{(∂_t)u[t,x]==(∂_(x,x))u[t,x],u[0,x]==0,u[t,0]==Sin[t],u[t,5]==0},u,{t,0,10},{x,0,5}]

Solve the Poisson equation over a (RefLink[Disk, paclet:ref/Disk]):

NDSolve[{-Laplacian[u[x,y],{x,y}]==1,DirichletCondition[u[x,y]==0,True]},u,{x,y}∈Disk[]]

Plot3D[Evaluate[u[x,y]/.%],{x,y}∈Disk[]]

Find a minimal surface over a (RefLink[Disk, paclet:ref/Disk]) with a sinusoidal boundary condition.

NDSolve[{-(InactiveDiv[((1/(Sqrt[1+(Subscript[∇, {x,y}])u[x,y].(Subscript[∇, {x,y}])u[x,y]]))) (InactiveGrad[u[x,y], {x,y}, Inactive]), {x,y}, Inactive])==0,DirichletCondition[u[x,y]==Sin[2π*(x+y)],True]},u,{x,y}∈Disk[]]

Plot3D[Evaluate[u[x,y]/.%],{x,y}∈Disk[]]

Solve a coupled nonlinear sine-Gordon equation over a region.

Ω=Polygon[2{{-Pi,-E},{0,-1},{0,0},{Pi,0},{Pi,E},{-Pi,E}}];
NDSolve[{Laplacian[u[x,y],{x,y}]+λ Cos[v[x,y]]  Sin[u[x,y]]==0 ,Laplacian[v[x,y],{x,y}]+λ Sin[u[x,y] ] Sin[v[x,y]]==0,DirichletCondition[{u[x,y]==1/10,v[x,y]==1/20},True]}/.λ->5,{u,v},{x,y}∈Ω]

Row[{Plot3D[Evaluate[u[x,y]/.%],{x,y}∈Ω,ImageSize->Medium],Plot3D[Evaluate[v[x,y]/.%],{x,y}∈Ω,ImageSize->Medium]}]

OpenersExampleSection1(Spacer1)Scope  (25)

OpenersExampleSubsection0(Spacer1)Ordinary Differential Equations  (8)

Specify any order equation. Reduction to normal form is done automatically:

NDSolve[{x''[t]+1/10x'[t]+Sin[x[t]]==1/2 Cos[t],x[0]==x'[0]==0},x,{t,0,100}]

Directly differentiate the solution to make a phase plot:

ParametricPlot[Evaluate[{x[t],x'[t]}/.%],{t,0,100},ColorFunction->Hue]

Directly specify a system of equations:

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,5}]

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

Solve for a vector-valued function:

s=NDSolve[{y''[x]+(({{1, 1, 1, 1}, {1, 2, 1, 2}, {1, 1, 3, 1}, {1, 2, 1, 4}})).y[x]==0,y[0]==y'[0]=={1,1,1,1}},y,{x,0,8}]

y[5]/.s

Plot the four components of the solution:

Plot[y[x] /. First[s],{x,0,8}]

Different equivalent ways of specifying a harmonic oscillator as a second-order equation:

NDSolve[{x''[t]+x[t]==0,x[0]==1,x'[0]==0},x,{t,10}]

Plot[Evaluate[{x[t],x'[t]}/.First[%]],{t,0,10}]

As a system of first-order equations:

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

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

Using a vector variable with the dimension deduced from the initial condition:

NDSolve[{z'[t]=={{0,1},{-1,0}}.z[t],z[0]=={1,0}},z,{t,0,10}]

Plot[z[t]/.First[%], {t,0,10}]

Use matrix-valued variables to compute the fundamental matrix solution:

A=RandomReal[{0,1},{5,5}];

fs=x/.First[NDSolve[{x'[t]==A.x[t],x[0]==IdentityMatrix[5]},x,{t,0,1}]]

Compare to the exact solution:

Plot[Norm[Flatten[fs[t]-MatrixExp[A t]]],{t,0,1}]

Define a Van der Pol equation:

vdp={x'[t]==y[t],y'[t]==-x[t]+1000(1-x[t]^2)y[t],x[0]==2,y[0]==0};

The solution's "stiff" behavior that the default solver automatically handles:

sol=NDSolve[vdp,{x,y},{t,2000}]

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

Other methods may not be able to solve this system:

NDSolve[{x'[t]==y[t],y'[t]==-x[t]+1000 (1-((x[t])^2)) y[t],x[0]==2,y[0]==0},{x,y},{t,2000},Method->"ExplicitRungeKutta"]

Equations can have multiple distinct solutions:

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

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

The solution y[x] is continuous, as it integrates the piecewise function once:

NDSolve[{y'[x]+Cos[y[x]]==Floor[x],y[0]==1},y,{x,0,3}]

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

The solution y[x] is differentiable, whereas y'[x] is continuous only:

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

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

OpenersExampleSubsection1(Spacer1)Partial Differential Equations  (5)

Nonlinear advection-diffusion equation in one dimension:

sol = NDSolve[{D[u[t,x], t] == 0.5 D[u[t,x],x,x]+u[t,x]D[u[t,x],x],  u[t,-Pi] == u[t, Pi] == 0 , u[0,x] == Sin[x]}, u, {t,0,2},{x,-Pi, Pi}]

Plot3D[Evaluate[u[t,x]/. First[sol]],{t,0,2},{x,-Pi, Pi}, PlotRange->All]

Define a system of PDEs of mixed parabolic-hyperbolic type:

pde={(∂_t)u[t,x]==(∂_x)((v[t,x]-1) (∂_x)u[t,x])+(16 x t-2 t-16 (v[t,x]-1)) (u[t,x]-1)+10 x (E^(-4 x)),(∂_t)v[t,x]==(∂_({x,2}))v[t,x]+(∂_x)u[t,x]+4 u[t,x]-4+(x^2)-2 t-10 t (E^(-4 x))};

bc={u[0,x]==1,v[0,x]==1,u[t,0]==1,v[t,0]==1,3 u[t,1]+Derivative[0,1][u][t,1]==3,5 Derivative[0,1][v][t,1]==(E^4) (u[t,1]-1)};

NDSolve[{pde,bc},{u,v},{x,0,1},{t,0,2}]

Plot3D[Evaluate[{u[t,x],v[t,x]}/.%],{x,0,1},{t,0,2},PlotStyle->{Directive[Opacity[0.7],Red],Directive[Opacity[0.7],Blue]}]

Nonlinear sine-Gordon equation in two spatial dimensions with periodic boundary conditions:

L = 4;
sol = 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[t,-L, y] == u[t, L, y], u[t, x, -L] == u[t, x, L], u[0,x, y] == Exp[-(x^2 + y^2)], Derivative[1,0,0][u][0,x,y] == 0}, u, {t,0,L/2}, {x,-L,L},{y,-L,L}]

Plot the solution at the final time:

Plot3D[u[L/2, x, y] /. First[sol], {x,-L,L},{y,-L,L}]

Plot the time evolution of a radial cross section of the solution:

Plot3D[u[t, x, 0] /. First[sol],{t,0, L/2}, {x,0,L}]

Solve a wave equation over a region with a slit:

Ω=RegionDifference[RegionDifference[Rectangle[{0,0},{2,1}],Rectangle[{9/10,0},{11/10,4/10}]],Rectangle[{9/10,6/10},{11/10,1}]];
sol=u/.First[NDSolve[{D[u[t,x,y],{t,2}]-Laplacian[u[t,x,y],{x,y}]==0,DirichletCondition[u[t,x,y]==0,True],u[0,x,y]==2*Exp[-125((x-0.25)^2+(y-0.5)^2)],Derivative[1,0,0][u][0,x,y]==0},u,{t,0,2},{x,y}∈Ω]]

ListAnimate[Table[Rasterize[Plot3D[sol[t,x,y],{x,y}∈Ω,PlotRange->{-0.75,2},AspectRatio->Automatic]],{t,0,2,1/25}],SaveDefinitions->True]

Solve a Poisson equation with periodic boundary conditions on curved boundaries:

Ω=RegionDifference[RegionUnion[Disk[],Rectangle[{0,-1},{2,1}]],Disk[{2,0}]];

ufun=First[u/.NDSolve[{-(Subscript[∇, {x,y}]^2)u[x,y]==1.,PeriodicBoundaryCondition[u[x,y],(x-2)^2+y^2==1,Function[x,x-{2,0}]],DirichletCondition[u[x,y]==0,(0<x<(2-10^-6)&&(y≤-1||y≥1))]},u,{x,y}∈Ω]]

Visualize the solution:

ContourPlot[ufun[x,y],{x,y}∈Ω,ColorFunction->"TemperatureMap",AspectRatio->Automatic]

OpenersExampleSubsection2(Spacer1)Boundary Value Problems  (5)

A nonlinear multipoint boundary value problem:

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

Plot[Evaluate[{x[t], y[t]} /. First[sol]],{t,0,4},Epilog->{Point[{0,1}],Point[{4,1}],Red,Point[{1,1}]}]

Solve a nonlinear diffusion equation ∇·(-(u^2) ∇u)=4 with Dirichlet and Neumann boundary conditions starting from an initial seed of u=1.

sol=NDSolve[{(InactiveDiv[(-((u[x])^2) (InactiveGrad[u[x], {x}, Inactive])), {x}, Inactive])==4+NeumannValue[2.,x==1],DirichletCondition[u[x]==1.,x==0]},u,{x}∈Line[{{0},{1}}],InitialSeeding->{u[x]==1}]

Visualize the result.

Plot[Evaluate[u[x] /. First[sol]],{x,0,1}]

Solve a nonlinear equation (∇^2)u+4u∇u=2 with Dirichlet boundary conditions starting from an initial seed of u=1.

sol=NDSolve[{(Subscript[∇, {x}]^2)u[x]+4u[x]D[u[x],x]==2,DirichletCondition[u[x]==1.,True]},u,{x}∈Line[{{0},{1}}],InitialSeeding->{u[x]==1}]

Visualize the result.

Plot[Evaluate[u[x] /. First[sol]],{x,0,1}]

Solve a complex-valued nonlinear reaction equation -(∇^2)u+i (u^2)=10 with Dirichlet boundary conditions.

sol=NDSolve[{-(Subscript[∇, {x}]^2)u[x]+I u[x]^2==10,DirichletCondition[u[x]==-6 I,x==0],DirichletCondition[u[x]==-Sqrt[2] I,x==1]},u,{x}∈Line[{{0},{1}}]]

Visualize the result.

ReImPlot[Evaluate[u[x] /. First[sol]],{x,0,1}]

Solve a boundary value problem with a nonlinear load term (∇^2)u=sin(u).

sol=NDSolve[{(Subscript[∇, {x}]^2)u[x]==Sin[u[x]],DirichletCondition[u[x]==5.,x==0],DirichletCondition[u[x]==5.1,x==π/2]},u,{x}∈Line[{{0},{π/2}}]]

Visualize the result.

Plot[Evaluate[u[x] /. First[sol]],{x,0,π/2}]

OpenersExampleSubsection3(Spacer1)Delay Differential Equations  (2)

Solve a delay differential with two constant delays and initial history function cos(t):

sol=NDSolve[{x'[t]==x[t](x[t-Pi]-x'[t-1]),x[t/;t≤0]==Cos[t]},x,{t,0,8}]

Discontinuities are propagated from t==0 at intervals equal to the delays:

Plot[Evaluate[{x[t],x'[t]} /. First[sol]],{t,0,8}, PlotRange->All]

Investigate stability for a linear delay differential equation:

Manipulate[
Module[{sol,y,t},
sol = First[y /. NDSolve[{y'[t] == λ y[t] + μ y[t - 1], y[t /; t ≤ 0]==1 - t}, y,
{t,0,10}]];
If[pp, ParametricPlot[{sol[t],sol[t-1]},{t,1,10},PlotRange->{{-3,3},{-3,3}}],
Plot[sol[t],{t,0,10}, PlotRange->{{0,10},{-3,3}}]]],{{pp,False, "Plot in phase plane"},{True, False}},{{λ, -1},-5,5},{{μ,1},-5,5}, FrameLabel->ToString[y'[t] == λ y[t]+μ y[t-1], TraditionalForm]]

OpenersExampleSubsection4(Spacer1)Hybrid and Discontinuous Equations  (4)

A differential equation with a discontinuous right-hand side using automatic event generation:

sol=NDSolve[{y'[t]==-Sign[y[t]],y[0]==1},y,{t,0,2}];

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

A differential equation whose right-hand side changes at regular time intervals:

sol=NDSolve[{y'[t]==a[t],y[0]==0,a[0]==1,WhenEvent[Mod[t,1]==0,a[t]->-a[t]]},y,{t,0,5},DiscreteVariables->a[t]];

Plot[y[t]/.sol,{t,0,5}]

Reflect a solution across the y axis each time it crosses the negative x axis:

de={x'[t]==y[t],y'[t]==-x[t]+.2 y[t]+1};
ic={x[0]==.5,y[0]==.5};

sol=NDSolve[{de,ic,WhenEvent[And[y[t]==0,x[t]<0],x[t]->-x[t]]},{x,y},{t,0,1000}];

ParametricPlot[{x[t],y[t]}/.sol,{t,0,200},PlotPoints->200]

Periodic solution with sliding mode:

Manipulate[Module[{x0, y0, sol},
{x0, y0} = p; sol =  First[NDSolve[{x'[t]==y[t],y'[t]==If[y[t] - 1>0, -1,.1 y[t] - x[t]],x[0]==x0,y[0]==y0},{x,y},{t,0,20}]];
ParametricPlot[{x[t], y[t]} /. sol,{t,0,20}, PlotRange->{{-1.5,1.5},{-1.5,1.5}}]],
{{p, {-1.5,1.5}},Locator}]

OpenersExampleSubsection5(Spacer1)Differential-Algebraic Equations  (1)

A differential equation with an algebraic constraint:

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

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

OpenersExampleSection2(Spacer1)Generalizations & Extensions  (1)

The names of functions need not be symbols:

eqns={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}]}

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

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

OpenersExampleSection3(Spacer1)Options  (29)

OpenersExampleSubsection6(Spacer1)AccuracyGoal and PrecisionGoal  (1)

Use defaults to solve a celestial mechanics equation with sensitive dependence on initial conditions:

eqn=z''[t]==-z[t]/(((((z[t])^2)+((((1+ Sin[2 π t]/2)/2))^2)))^(3/2));

Plot[Evaluate[z[t]/.NDSolve[{eqn,z[0]==1,z'[0]==0},z,{t,0,40}]],{t,0,40}]

Higher accuracy and precision goals give a different result:

Plot[Evaluate[z[t]/.NDSolve[{eqn,z[0]==1,z'[0]==0},z,{t,0,40},AccuracyGoal->10,PrecisionGoal->10]],{t,0,40}]

Increasing the goals extends the correct solution further:

Plot[Evaluate[z[t]/.NDSolve[{eqn,z[0]==1,z'[0]==0},z,{t,0,40},AccuracyGoal->20,PrecisionGoal->20,WorkingPrecision->35]],{t,0,40}]

OpenersExampleSubsection7(Spacer1)DependentVariables  (1)

Set up a very large system of equations:

n = 1000;
vars = Table[(Subscript[x, i])[t], {i, n}];
eqns = Table[j = Mod[i, n] + 1;{(Subscript[x, i])'[t] == 1/((((Subscript[x, i])[t] + (Subscript[x, j])[t]))^2), (Subscript[x, i])[0] == 1/i}, {i, n}];

Short[eqns, 3]

Solve for all the dependent variables, but save only the solution for (Subscript[x, 1]):

sol = NDSolve[eqns, (Subscript[x, 1]), {t,0,100}, DependentVariables->vars]

Plot[(Subscript[x, 1])[t] /. sol, {t,0,100}, PlotRange->All]

OpenersExampleSubsection8(Spacer1)EvaluationMonitor  (2)

Total number of evaluations:

Module[{c=0},NDSolve[{y'[x]==y[x]Cos[x+y[x]],y[0]==1},y,{x,0,30},EvaluationMonitor:>c++];c]

The distance between successive evaluations; negative distance means a rejected step:

ListLinePlot[Differences[Reap[NDSolve[{y'[x]==y[x]Cos[x+y[x]],y[0]==1},y,{x,0,30},EvaluationMonitor:>Sow[x]]][[2,1]]]]

OpenersExampleSubsection9(Spacer1)InitialSeedings  (2)

Specify an initial seeding of 0 for a boundary value problem:

NDSolveValue[{-2 u[x] u'[x]-u''[x]==1,DirichletCondition[u[x]==0,x==0],DirichletCondition[u[x]==0,x==1]},u,{x}∈Line[{{0},{1}}],InitialSeeding->{u[x]==0}]

Specify an initial seeding that depends on a spatial coordinate:

NDSolveValue[{-2 u[x] u'[x]-u''[x]==0,DirichletCondition[u[x]==0,x==0],DirichletCondition[u[x]==1,x==1]},u,{x}∈Line[{{0},{1}}],InitialSeeding->{u[x]==x}]

OpenersExampleSubsection10(Spacer1)InterpolationOrder  (1)

Use (RefLink[InterpolationOrder, paclet:ref/InterpolationOrder])->(RefLink[All, paclet:ref/All]) to get interpolation the same order as the method:

Timing[a = NDSolve[{x'[t] == y[t], y'[t]== -x[t],x[0]==1, y[0]==0},{x,y},{t,0,1000},InterpolationOrder->All]]

This is more time-consuming than the default interpolation order used:

Timing[d = NDSolve[{x'[t] == y[t], y'[t]== -x[t],x[0]==1, y[0]==0},{x,y},{t,0,1000}]]

It is much better in between steps:

Plot[Evaluate[{First[x[t] /. a],First[x[t]/. d]} - Cos[t]],{t,0,10}, PlotRange->All, PlotStyle->{Directive[Thick, Blue],Red}]

OpenersExampleSubsection11(Spacer1)MaxStepFraction  (1)

Features with small relative size in the integration interval can be missed:

Table[Plot[Evaluate[y[x]/.NDSolve[{y''[x]+y[x] ==Exp[-(1000 (x/L - 1/3))^2], y[0]==y'[0]==0},y,{x,L}]], {x,L/4,L/2},PlotRange->All],{L, {1, 10, 100, 1000}}]

Use (RefLink[MaxStepFraction, paclet:ref/MaxStepFraction]) to ensure features are not missed, independent of interval size:

Table[Plot[Evaluate[y[x]/.NDSolve[{y''[x]+y[x]==Exp[-(1000 (x/L - 1/3))^2], y[0]==y'[0]==0},y,{x,L},MaxStepFraction->0.001]], {x,L/4,L/2},PlotRange->All],{L, {1, 10, 100, 1000}}]

OpenersExampleSubsection12(Spacer1)MaxSteps  (1)

Integration stops short of the requested interval:

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

More steps are needed to resolve the solution:

sol = NDSolve[{y''[x] + x y[x]  == 0, y[0] == 1, y'[0] == 0}, y, {x,0,200}, MaxSteps->10^6]

Plot the solution in the phase plane:

ParametricPlot[Evaluate[{y[x], y'[x]} /. First[sol]], {x, 0,200}, ColorFunction->Function[Hue[#3]], AspectRatio->1]

For an infinite integration of an oscillator, a maximum number of steps is reached:

NDSolve[{x''[t]+x[t]==0,x[0]==1,x'[0]==0},x,{t,Infinity}]

More steps can be requested:

NDSolve[{x''[t]+x[t]==0,x[0]==1,x'[0]==0},x,{t,Infinity},MaxSteps->20000]

OpenersExampleSubsection13(Spacer1)MaxStepSize  (2)

The default step control may miss a suddenly varying feature:

Plot[Evaluate[(y[x] /. NDSolve[{y''[x] + (1 +  Sech[1000 (x - Pi)]) y[x] == 0, y[0] == 1, y'[0] == 0}, y,{x,0,10}]) - Cos[x]], {x,0,10}]

A smaller (RefLink[MaxStepSize, paclet:ref/MaxStepSize]) setting ensures that (RefLink[NDSolve, paclet:ref/NDSolve]) catches the feature:

Plot[Evaluate[(y[x] /. NDSolve[{y''[x] + (1 +  Sech[1000 (x - Pi)]) y[x] == 0, y[0] == 1, y'[0] == 0}, y,{x,0,10}, MaxStepSize->0.01]) - Cos[x]], {x,0,10}]

Attempting to compute the number of positive integers less than (e^5) misses several events:

Block[{n = 0},NDSolve[{y'[t] == y[t], y[-1] == (E^(-1)),WhenEvent[Sin[π y[t]]==0, n++]}, y, {t,5}];n]

Setting a small enough (RefLink[MaxStepSize, paclet:ref/MaxStepSize]) ensures that none of the events are missed:

Block[{n = 0},NDSolve[{y'[t] == y[t], y[-1] == (E^(-1)),WhenEvent[Sin[π y[t]]==0, n++]}, y, {t,5},MaxStepSize->0.001];n]

OpenersExampleSubsection14(Spacer1)Method  (12)

OpenersExampleSubsubsection0(Spacer1)TimeIntegration  (4)

Specify an explicit Runge–Kutta method to be used for the time integration of a differential equation:

NDSolve[{y'[x]==y[x]Cos[x+y[x]],y[0]==1},y,{x,0,30},Method->{"TimeIntegration"->"ExplicitRungeKutta"}]

Specify an explicit Runge–Kutta method of order 8 to be used for the time integration:

NDSolve[{y'[x]==y[x]Cos[x+y[x]],y[0]==1},y,{x,0,30},Method->{"TimeIntegration"->{"ExplicitRungeKutta","DifferenceOrder"->8}}]

Specify an explicit Euler method to be used for the time integration of a differential equation:

NDSolve[{y''[t]==-y[t],y[0]==1,y'[0]==0},y,{t,0,1},Method->{"TimeIntegration"->"ExplicitEuler"}]

Extrapolation tends to take very large steps:

ListLinePlot[Differences[Reap[NDSolve[{y'[x]==y[x]Cos[x+y[x]],y[0]==1},y,{x,0,30},StepMonitor:>Sow[x],Method->{"TimeIntegration"->"Extrapolation"}]][[2,1]]]]

OpenersExampleSubsubsection1(Spacer1)PDEDiscretization  (2)

Solutions of Burgers' equation may steepen, leading to numerical instability:

s1=NDSolve[{(∂_t)u[t,x]==(1/1000)(∂_(x,x))u[t,x]- u[t,x](∂_x)u[t,x],u[0,x]==Sin[2 π x],u[t,0]==u[t,1]},u,{t,0,2},{x,0,1}]

end = Part[u  /. s1, 1, 1, 1,-1];
Plot[u[end,x] /. s1, {x,0,1}]

Specify a spatial discretization sufficiently fine to resolve the front:

s2=NDSolve[{(∂_t)u[t,x]==(1/1000)(∂_(x,x))u[t,x]- u[t,x](∂_x)u[t,x],u[0,x]==Sin[2 π x],u[t,0]==u[t,1]},u,{t,0,2},{x,0,1}, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->1000}}}]

end = Part[u  /. s2, 1, 1, 1,-1];
Plot[u[end,x] /. s2, {x,0,1}]

After the front forms, the solution decays relatively rapidly:

X = Part[u /. s2, 1, 3, 2];
Plot[Max[u[t, X] /. s2], {t,0,2}]

Specify use of the finite element method for spatial discretization:

s2=NDSolve[{(∂_t)u[t,x]==(1/1000)(∂_(x,x))u[t,x],u[0,x]==Sin[2 π x],DirichletCondition[u[t,x]==0,True]},u,{t,0,2},{x,0,1}, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"FiniteElement"}}}]

OpenersExampleSubsubsection2(Spacer1)BoundaryValues  (1)

Solve a boundary value problem:

sol0 = First[NDSolve[{x''[t] + Sin[x[t]] == 0, x[0] == x[10] == 0},x,t]]

With the default option, the method finds the trivial solution:

Table[x[t] /. sol0, {t,0,10}]

Specify different starting conditions for the "Shooting" method to find different solutions:

sols = Map[First[NDSolve[{x''[t] + Sin[x[t]] == 0, x[0] == x[10] == 0},x,t, Method->"BoundaryValues"->{"Shooting", "StartingInitialConditions"->{x[0] == 0, x'[0] == #}}]]&,
{1.5, 1.75, 2}];
Plot[Evaluate[x[t] /. sols], {t,0,10}, PlotStyle->{Black, Blue, Green}]

OpenersExampleSubsubsection3(Spacer1)DiscontinuityProcessing  (1)

(RefLink[NDSolve, paclet:ref/NDSolve]) automatically does processing for discontinuous functions like (RefLink[Sign, paclet:ref/Sign]):

sol1=NDSolve[{x'[t]==Sign[1-x[t]],x[0]==0},x,{t,0,2}]

If the processing is turned off, (RefLink[NDSolve, paclet:ref/NDSolve]) may fail at the discontinuity point:

NDSolve[{x'[t]==Sign[1-x[t]],x[0]==0},x,{t,0,2}, Method->{"DiscontinuityProcessing"->False}]

With some time integration methods, the solution may be very inaccurate:

sol2 = NDSolve[{x'[t]==Sign[1-x[t]],x[0]==0},x,{t,0,2}, Method->{"DiscontinuityProcessing"->False, "TimeIntegration"->"Extrapolation"}]

Plot[{x[t]/.sol1 ,x[t]/.sol2}, {t, 0,2}]

An equivalent way to find the solution is to use "DiscontinuitySignature":

sol3=NDSolve[{x'[t]==v[t],WhenEvent[1==x[t],v[t]->"DiscontinuitySignature"],x[0]==0,v[0]==1},{x,v},{t,0,2}, DiscreteVariables->{Element[v,{-1,0,1}]}]

The solutions are effectively identical:

Table[(x[t]/.sol1) -(x[t]/.sol3), {t, 0,2,.25}]

The discontinuity signature is 0 when the solution is in sliding mode:

Plot[v[t] /. sol3,{t,0,2}]

OpenersExampleSubsubsection4(Spacer1)EquationSimplification  (2)

The solution cannot be completed because the square root function is not sufficiently smooth:

NDSolve[{x'[t]^2+x[t]^2==1,x[0]==1/2},x,{t,0,10 Pi}]

One solution can be found by forming a residual and solving as a DAE system:

s1 = NDSolve[{x'[t]^2+x[t]^2==1,x[0]==1/2},x,{t,0,10 Pi}, Method->{"EquationSimplification"->"Residual"}]

The other solution branch can be given by specifying a consistent value of x':

s2 = NDSolve[{x'[t]^2+x[t]^2==1,x[0]==1/2, x'[0]==-Sqrt[3/4]},x,{t,0,10 Pi}, Method->{"EquationSimplification"->"Residual"}]

Plot[Evaluate[x[t] /. {s1,s2}], {t, 0, 10 Pi}]

With the suboption "SimplifySystem"->(RefLink[True, paclet:ref/True]), (RefLink[NDSolve, paclet:ref/NDSolve]) uses symbolic solutions for components with a sufficiently simple form:

sol = NDSolve[{x'[t]+y[t]==Sin[t],y'[t]+z[t]==Sin[t],y[t]==Cos[t],x[0]==1},{x,y,z},{t,0,2Pi},Method->{"EquationSimplification"->{Automatic,
"SimplifySystem"->True}}]

Plot[Evaluate[{x[t],y[t],z[t]}/.sol],{t,0,2Pi},
PlotLegends->{x[t],y[t],z[t]}]

OpenersExampleSubsubsection5(Spacer1)IndexReduction  (1)

An index 3 formulation of a constrained pendulum using index reduction:

eqns={x''[t]==T[t]x[t],y''[t]==T[t] y[t]-1,x[t]^2+y[t]^2==1};
init={x[0]==1,y[0]==0,y'[0]==-14/10};
vars={x,y, T};

The default method can only solve index 1 problems:

NDSolve[{eqns, init}, vars, {t,0,10}]

The problem resulting from symbolic index reduction can be solved:

s1 = First[NDSolve[{eqns,init},vars,{t,0,10},Method->{"IndexReduction"->True}]]

Solve using reduction to index 0 and a projection method to maintain the constraints:

s2 = First[NDSolve[{eqns,init},vars,{t,0,10},Method->{"IndexReduction"->{True, "ConstraintMethod"->"Projection"}}]]

Plot implicit energy constraint for the two solutions at the time steps:

energy[sol_] := Module[{times = Part[x /. sol, 3, 1], es}, 
es[t_] = 1+y[t]+(1/2) (((-y[t] x'[t]+x[t] y'[t]))^2) /. sol;Transpose[{times,  es[times] - es[0]}]];

ListPlot[{energy[s1], energy[s2]}]

OpenersExampleSubsubsection6(Spacer1)DAEInitialization  (1)

Use forward collocation for initialization to avoid problems with the (RefLink[Abs, paclet:ref/Abs]) term at 0.

NDSolve[{x'[t]==Abs[u[t]],u[t]==Sin[t],y[t]==x[t],y[0]==10},{x,y,u},{t,0,10}, Method->{"DAEInitialization"->{"Collocation", "CollocationDirection"->"Forward"}}]

OpenersExampleSubsection15(Spacer1)NormFunction  (1)

Plot the actual solution error when using different error estimation norms:

L = 10;Table[sol[p] = 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] == -Exp[-x^2], u[t,-L] == u[t,L]}, u, {t,0,4 L},{x,-L,L}, NormFunction->(Norm[#, p]&),Method->{"MethodOfLines","SpatialDiscretization"->{"TensorProductGrid","MaxPoints"->125,"MinPoints"->125,"DifferenceOrder"->"Pseudospectral"}}];
Plot[Evaluate[(u[4L,x] /. sol[p]) -  (Exp[-x^2] - 2(Sqrt[π]))],{x,-L,L}, PlotRange->All],{p, {1,2,3,∞}}]

A plot of the best solution:

Plot3D[u[t,x] /. sol[1], {t,0,40},{x,-L,L}, Mesh->False, PlotPoints->20]

OpenersExampleSubsection16(Spacer1)StartingStepSize  (1)

For a very large interval, a short-lived feature near the start may be missed:

feature[t_?NumberQ] := If[t < 1, t Sin[10 Pi t], 0];sol = NDSolve[{x''[t] + 10^-4 x[t] == feature[t],x[0] == x'[0] == 0}, x,  {t, 0, 10000}]

Max[Abs[x[RandomReal[{0, 10000}, 10000]] /. sol]]

Setting a sufficiently small step size to start with ensures that the input is not missed:

sol = NDSolve[{x''[t] + 10^-4 x[t] == feature[t],x[0] == x'[0] == 0}, x,  {t, 0, 10000}, StartingStepSize->0.1]

{Plot[x[t] /. sol, {t,0,2}],Plot[x[t] /. sol, {t, 0, 10000}]}

OpenersExampleSubsection17(Spacer1)StepMonitor  (3)

Plot the solution at each point where a step is taken in the solution process:

ListLinePlot[Reap[NDSolve[{y'[x]==y[x]Cos[x+y[x]],y[0]==1},y,{x,0,30},StepMonitor:>Sow[{x,y[x]}]]][[2,1]]]

Total number of steps involved in finding the solution:

Module[{c=0},NDSolve[{y'[x]==y[x]Cos[x+y[x]],y[0]==1},y,{x,0,30},StepMonitor:>c++];c]

Differences between values of x at successive steps:

ListLinePlot[Differences[Reap[NDSolve[{y'[x]==y[x]Cos[x+y[x]],y[0]==1},y,{x,0,30},StepMonitor:>Sow[x]]][[2,1]]]]

OpenersExampleSubsection18(Spacer1)WorkingPrecision  (1)

Error in the solution to a harmonic oscillator over 100 periods:

T = 200 Pi;
Timing[y[T] /. NDSolve[{y''[x] + y[x] == 0, y[0] == 0, y'[0] == 1},y,{x,0,T}, MaxSteps->Infinity]]

When the working precision is increased, the local tolerances are correspondingly increased:

Timing[y[T] /. NDSolve[{y''[x] + y[x] == 0, y[0] == 0, y'[0] == 1},y,{x,0,T}, WorkingPrecision->32, MaxSteps->10^6]]

With a large working precision, sometimes the "Extrapolation" method is quite effective:

Timing[y[T] /. NDSolve[{y''[x] + y[x] == 0, y[0] == 0, y'[0] == 1},y,{x,0,T}, WorkingPrecision->32, Method->"Extrapolation"]]

Error in the solution to a harmonic oscillator over 100 periods:

T = 200 Pi;
Timing[y[T] /. NDSolve[{y''[x] + y[x] == 0, y[0] == 0, y'[0] == 1},y,{x,0,T}]]

When the working precision is increased, the local tolerances are correspondingly increased:

Timing[y[T] /. NDSolve[{y''[x] + y[x] == 0, y[0] == 0, y'[0] == 1},y,{x,0,T}, WorkingPrecision->32,MaxSteps->Infinity]]

With a large working precision, sometimes the "Extrapolation" method is quite effective:

Timing[y[T] /. NDSolve[{y''[x] + y[x] == 0, y[0] == 0, y'[0] == 1},y,{x,0,T}, WorkingPrecision->32, Method->"Extrapolation"]]

OpenersExampleSection4(Spacer1)Applications  (31)

OpenersExampleSubsection19(Spacer1)Ordinary Differential Equations  (5)

Simulate Duffing's equation for a particle in a double potential well:

s=NDSolve[{x''[t]+0.15x'[t]-x[t]+x[t]^3==0.3Cos[t],x[0]==-1,x'[0]==1},x,{t,0,50}]

Plot[Evaluate[x[t]/.s],{t,0,50}]

The solution depends strongly on initial conditions:

s=NDSolve[{x''[t]+0.15x'[t]-x[t]+x[t]^3==0.3Cos[t],
 x[0]==-1,x'[0]==1.001},x,{t,0,50}]

Plot[Evaluate[x[t]/.s],{t,0,50}]

The Lorenz equations [(WebLink[more info, http://mathworld.wolfram.com/LorenzEquations.html])]:

s=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}]

ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.s],{t,0,200},PlotPoints->1000,ColorFunction->(Hue[#4]&)]

The Lotka–Volterra predator-prey equations [(WebLink[more info, http://mathworld.wolfram.com/Lotka-VolterraEquations.html])]:

NDSolve[{y'[t]==y[t](x[t]-1),x'[t]==x[t](2-y[t]),x[0]==1,y[0]==2.7},{x,y},{t,0,10}]

Phase plane plot:

ParametricPlot[Evaluate[{x[t],y[t]}/.First[%]],{t,0,10}]

Look at the appearance of the blue sky catastrophe orbit in the Gavrilov–Shilnikov model:

eqns = {
x'[t]==x[t](2+μ-10(x[t]^2+y[t]^2)) + z[t]^2 + y[t]^2 +2 y[t],
y'[t]==-z[t]^3-(1+y[t])(z[t]^2+y[t]^2+2 y[t])-4 x[t]+μ y[t],
z'[t]==(1+y[t]) z[t]^2 +x[t]^2-ϵ};

Table[Block[{ϵ=0.0357,T=1000},
sol = First[NDSolve[{eqns, x[0]==y[0]==z[0]==1},{x,y,z},{t,0,T}]];
ParametricPlot3D[Evaluate[{x[t],y[t],z[t]} /. sol],{t,T/4,T}, BoxRatios->1,PlotPoints->1000,PlotRange->All]],{μ,.455,.456,.001}]

Reduced 3-body problem (WebLink[[more info], http://www.wolframscience.com/nksonline/page-973a-text]):

eqn=z''[t]==-((z[t])/(((((z[t])^2)+((((1/2) (1+0.1 Sin[2 Pi t])))^2)))^(3/2)));

NDSolve[{eqn,z[0]==1,z'[0]==0},z,{t,0,30}]

Plot[Evaluate[z[t]/.%],{t,0,30}]

A formulation suitable for a number of different initial conditions:

eqn=D[z[t,s],t,t]==-((z[t,s])/(((((z[t,s])^2)+((((1/2) (1+0.1 Sin[2 Pi t])))^2)))^(3/2)));

sol = NDSolve[{eqn,z[0,s]==s,Derivative[1,0][z][0,s]==0},z,{t,0,10}, {s,1,2}]

Plot3D[Evaluate[z[tp, z0]/.sol],{tp,0,10},{z0,1,2}]

OpenersExampleSubsection20(Spacer1)Partial Differential Equations  (6)

Simple model for soil temperature at depth x with periodic heating at the surface:

NDSolve[{D[u[t,x],t]==D[u[t,x],x,x],u[0,x]==0,u[t,0]==Sin[t],u[t,5]==0},u,{t,0,10},{x,0,5}]

Plot3D[Evaluate[u[t,x]/.%],{t,0,10},{x,0,5}, PlotRange->All]

Simple wave evolution with periodic boundary conditions:

s=NDSolve[{(∂_(t,t))u[t,x]==(∂_(x,x))u[t,x],u[0,x]==(E^(-(x^2))),Derivative[1,0][u][0,x]==0,u[t,-10]==u[t,10]},u,{t,0,40},{x,-10,10}]

Plot the solution:

DensityPlot[Evaluate[First[u[t,x]/.s]],{t,0,40},{x,-10,10},PlotPoints->30]

Wolfram's nonlinear wave equation [(WebLink[more info, http://www.wolframscience.com/nksonline/page-923a-text])]:

s=NDSolve[{(∂_(t,t))u[t,x]==(∂_(x,x))u[t,x]+(1-((u[t,x])^2)) (1+2 u[t,x]),u[0,x]==(E^(-(x^2))),Derivative[1,0][u][0,x]==0,u[t,-10]==u[t,10]},u,{t,0,10},{x,-10,10}]

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

DensityPlot[Evaluate[u[-t,x]/.s],{x,-10,10},{t,0,-10}]

Wolfram's nonlinear wave equation in two space dimensions:

NDSolve[{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]),u[0,x,y]==(E^(-((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}]

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

A soliton profile perturbed by a periodic potential in a nonlinear Schrödinger equation:

L=50;
sol=NDSolve[{I D[u[t,x],t]+D[u[t,x],x,x]+2 Abs[u[t,x]]^2 u[t,x]==0.1 (1-Cos[Pi x/L]),u[0,x]==Sech[x] Exp[I x],u[t,-L]==u[t,L]},u,{t,0,10},{x,-L,L}]

Plot3D[Evaluate[Abs[u[t,x]/.First[sol]]],{t,0,10},{x,-L,L},PlotPoints->60,MaxRecursion->3,Mesh->None]

Use Stokes's equation to compute the fluid velocity field in a narrowing channel:

Ω=RegionUnion[Rectangle[{0,0},{1,1/2}],Rectangle[{1,1/10},{2,2/5}]];
{xVel,yVel}={u,v}/.First[NDSolve[{{(InactiveDiv[({{-μ,0},{0,-μ}}.(InactiveGrad[u[x,y], {x,y}, Inactive])), {x,y}, Inactive])+Derivative[1,0][p][x,y],(InactiveDiv[({{-μ,0},{0,-μ}}.(InactiveGrad[v[x,y], {x,y}, Inactive])), {x,y}, Inactive])+Derivative[0,1][p][x,y],Derivative[1,0][u][x,y]+Derivative[0,1][v][x,y]}=={0,0,0}/.μ->1,{
DirichletCondition[u[x,y]==4*0.3*y*(0.5-y)/(0.41)^2,x==0.],DirichletCondition[{u[x,y]==0.,v[x,y]==0.},0<x<2],
DirichletCondition[p[x,y]==0.,x== 2]}},{u,v,p},{x,y}∈Ω,Method->{"FiniteElement","InterpolationOrder"->{u->2,v->2,p->1},"MeshOptions"->{"MaxCellMeasure"->0.0005}}]];

VectorPlot[{xVel[x,y], yVel[x,y]},{x,y}∈Ω, AspectRatio->Automatic,StreamPoints->6,StreamColorFunctionScaling->False]

OpenersExampleSubsection21(Spacer1)Delay Differential Equations  (1)

View solutions of the Mackey–Glass delay differential equation for respiratory dynamics:

Manipulate[
Module[{sol,x,t},
sol = First[ NDSolve[{x'[t] == a x[t - τ]/(1 + x[t - τ]^10) -b  x[t], x[t /; t ≤ 0]==1/2},x,{t,0,500}]];
ParametricPlot[Evaluate[{x[t], x[t - τ]} /. sol],{t,300,500}, PlotRange->{{0,2},{0,2}}]],{{a,.2},0,1},{{b,.1},0,1},
{{τ,17},1,20}]

OpenersExampleSubsection22(Spacer1)Hybrid Differential Equations  (5)

Simulate a bouncing ball that retains 95% of its velocity in each bounce:

NDSolve[{y''[t]==-9.81,y[0]==5,y'[0]==0,WhenEvent[y[t]==0,y'[t]->-0.95y'[t]]},y,{t,0,10}];

Plot[y[t]/.%,{t,0,10}]

Model a ball bouncing down steps:

ballsteps=NDSolve[{x''[t]==0,y''[t]==-9.8,y[0]==6,y'[0]==0,x[0]==0,x'[0]==1,a[0]==5,WhenEvent[Mod[x[t],1]==0,If[a[t]>0,a[t]->a[t]-1,"RemoveEvent"]],WhenEvent[y[t]==a[t],{x'[t],y'[t]}->.9{x'[t], -y'[t]}]},{x,y,a},{t,0,15},DiscreteVariables->{a}];

Show[ParametricPlot[Evaluate[{{x[t],y[t]},{x[t],a[t]}} /. ballsteps],{t,0,15}], Plot[{0,Floor[6-x]},{x,-1,15}, Filling->{2->0}, Exclusions->None],Frame->{{True,False},{True,False}},FrameLabel->{"x","y"}]

Each time a linear oscillator solution crosses the negative x axis, reflect it across the y axis:

de[δ_]:={(Subscript[x, 1])'[t]==(Subscript[x, 2])[t],(Subscript[x, 2])'[t]==-(Subscript[x, 1])[t]+2δ (Subscript[x, 2])[t]+1};
ic={(Subscript[x, 1])[0]==.5,(Subscript[x, 2])[0]==.5};

{sol,{points}}=NDSolve[{de[0.1],ic,WhenEvent[And[(Subscript[x, 2])[t]==0,(Subscript[x, 1])[t]<0],{(Subscript[x, 1])[t]->-(Subscript[x, 1])[t],Sow[-(Subscript[x, 1])[t]]}]},{(Subscript[x, 1]),(Subscript[x, 2])},{t,0,1000}]//Reap;

The solution of this reset oscillator exhibits chaotic behavior:

Map[ParametricPlot[{(Subscript[x, 1])[t], (Subscript[x, 2])[t]}/.sol, {t, 0, #}, PlotPoints->500,PlotRange->{{-1,3},{-2,2}}]&,{10,100,400}]

Plot the solution on the negative x axis with a histogram of the reflection points:

Grid[{{Histogram[points,{-1,0, .05}]},{ParametricPlot[{(Subscript[x, 1])[t], (Subscript[x, 2])[t]}/.sol, {t, 0, 1000}, PlotRange->{{-1,0},{-0.3,0}},PlotPoints->2000,Axes->None]}}]

Model a one-degree-of-freedom impact oscillator with sinusoidal forcing:

sol[ω_]:=With[{r=0.95},
NDSolve[{u'[t]==v[t],v'[t]==-u[t]+Cos[ω t],u[0]==0,v[0]==0,WhenEvent[u[t]==0,v[t]->-r v[t]]},{u, v},{t,0,300}]];

Table[ParametricPlot[Evaluate[{u[t],v[t]}/.sol[ω]],{t,200,300}, PlotPoints->200,PlotRange->{{-.1,1},{-1,1}}],{ω,{3,2.76,2.9}}]

Model a damped oscillator that gets a kick at regular time intervals:

system={x''[t]+.1x'[t]+x[t]==0,x[0]==1,x'[0]==0};
control=WhenEvent[Mod[t,1],x'[t]->x'[t]+1];

sol=NDSolve[{system,control},x,{t,0,100}]

Plot[Evaluate[{x[t],x'[t]}/.sol],{t,0,100},PlotPoints->120]

The trajectory eventually settles into a consistent orbit:

ParametricPlot[{x[t],x'[t]}/.sol,{t,90,100},AspectRatio->1]

OpenersExampleSubsection23(Spacer1)Mechanical Systems  (3)

Model the motion of a pendulum in Cartesian coordinates. Derive the governing equations using Newton's second law of motion, m (x')'(t)==∑(Subscript[F, x])  and m (y')'(t)==∑(Subscript[F, y]), with a force diagram:

deqns={x''[t]==λ[t] x[t],y''[t]==λ[t] y[t]-9.81};
aeqns={((x[t])^2)+((y[t])^2)==(1^2)};
ics={x[0]==1,y'[0]==1};

Simulate the system:

sol1=NDSolve[{deqns,aeqns,ics},{x,y,λ},{t,0,5},Method -> {"IndexReduction" -> True}];

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

Add damping to the pendulum so that it slows down over time:

deqns2={x''[t]==λ[t]x[t]-0.5 x'[t], y''[t]==λ[t] y[t]-9.8-0.5 y'[t]};

sol2=NDSolve[{deqns2,aeqns,ics},{x,y,λ},{t,0,20},Method -> {"IndexReduction" -> True}];

The pendulum stabilizes at the vertical fixed point:

Plot[Evaluate[{x[t],y[t]}/.sol2],{t,0,20}]

Change the rod to a stiff spring by modifying the constraint ((x(t))^2)+((y(t))^2)==(l^2):

aeqns3={Sqrt[((x[t])^2)+((y[t])^2)](1+λ[t] κ)==1};
κ=5 10^-4;

sol3=Quiet@NDSolve[{deqns,aeqns3,ics},{x,y,λ},{t,0,5},Method->{"IndexReduction"->True},AccuracyGoal->6];

The solution contains high-frequency spring oscillations:

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

Model a double pendulum of unit mass and length that is released from the horizontal plane. Begin by deriving the equations of motion using Newton's second law:

deqns={x1''[t]==λ1[t] x1[t]-λ2[t](x2[t]-x1[t]), y1''[t]==λ1[t] y1[t]-λ2[t](y2[t]-y1[t])-9.81,
x2''[t]==λ2[t](x2[t]-x1[t]),y2''[t]==λ2[t] (y2[t]-y1[t])- 9.81};
aeqns={((x1[t])^2)+((y1[t])^2)==1,(((x2[t]-x1[t]))^2)+(((y2[t]-y1[t]))^2) ==1};
ics={x1[0]==1,y1'[0]==1,x2[0]==2,y2'[0]==1};

Simulate the system by enforcing the equations and constraints as invariants:

sol=Quiet@NDSolve[{deqns,aeqns,ics},{x1,y1,x2,y2,λ1,λ2},{t,0,15},Method->{"IndexReduction"->{True,"ConstraintMethod"->"Projection"}}] ;

Visualize the motion of the double pendulum:

Animate[Graphics[{{Thickness[0.02],Red,Line[{{0,0},{x1[t],y1[t]}}]},{Thickness[0.02],Blue,Line[{{x1[t],y1[t]},{x2[t],y2[t]}}]}}/.sol,PlotRange->{{-2.1,2.1},{0.6,-2.1}},Frame -> True],{t,0,15},SaveDefinitions->True,AnimationRunning->False]

Model a block on a moving conveyor belt anchored to a wall by a spring using different models for the friction force (Subscript[F, f]), including viscous, Coulomb, Stribeck, and static. Compare positions and velocities for the different models:

m=l=1.;
beltv[t_]=.1;
spring[x_]=1000.(l-x);

Newton's equation for the block:

sys:={m x''[t]==spring[x[t]]+friction[x'[t]],x[0]==1, x'[0]==0};

Viscous friction is proportional to the relative velocity (Subscript[F, visc])==-α v:

viscous[v_]:=-30.(v-beltv[t]);
friction[v_]:=viscous[v];

{pos1,vel1}=NDSolveValue[sys,{x[t],x'[t]},{t,0,2}];

The block stabilizes just above the spring's natural length of 1:

{Plot[{pos1},{t,0,1},PlotRange->All],Plot[{vel1,beltv[t]},{t,0,1},PlotRange->All]}

Coulomb friction is proportional to the sign of relative velocity (Subscript[F, coul])==β sgn(v):

coulomb[v_]:=-25.Sign[v-beltv[t]];
friction[v_]:=viscous[v]+coulomb[v];

{pos2,vel2}=NDSolveValue[sys,{x[t],x'[t]},{t,0,20}];

The block moves with the belt until the spring force is strong enough:

{Plot[{pos2},{t,0,1},PlotRange->All],Plot[{vel2,beltv[t]},{t,0,1},PlotRange->All]}

Stribeck friction is a refined Coulomb friction (Subscript[F, str])==γ sgn(v) (e^(-2 (Abs[v]))):

stribeck[v_]:=-.3Sign[v]Exp[-2Abs[v]];
friction[v_]:=viscous[v]+coulomb[v]+stribeck[v];

{pos3,vel3}=NDSolveValue[sys,{x[t],x'[t]},{t,0,20}];

The variation at low velocities is slightly reduced:

{Plot[{pos3},{t,0,1},PlotRange->All],Plot[{vel3,beltv[t]},{t,0,1},PlotRange->All]}

Static friction holds the block in place until the spring force exceeds some value μ depending on roughness of surfaces. Use the discrete variable stuck set to 1 when the block is stuck and 0 otherwise:

friction[v_]:=viscous[v]+coulomb[v]+stribeck[v];

staticsys:={x''[t]==If[stuck[t]==1,beltv'[t],spring[x[t]]+friction[x'[t]]],x[0]==1, x'[0]==0};

Check whether the spring force is smaller than μ, and if the block is not moving relative to the belt:

μ=100;
stick=WhenEvent[x'[t]==beltv[t],stuck[t]->Boole[((spring[x[t]])^2)<(μ^2)]];
slip=WhenEvent[((spring[x[t]])^2)>(μ^2),stuck[t]->0];

{pos4,vel4}=NDSolveValue[{staticsys,stick,slip,stuck[0]==0},{x[t],x'[t]},{t,0,8},DiscreteVariables->stuck[t]];

The block repeatedly sticks to the belt, then slips away due to the spring force:

{Plot[{pos4},{t,0,2},PlotRange->All],Plot[{vel4,beltv[t]},{t,0,2},PlotPoints->100,PlotRange->All]}

Compare the different models:

Plot[{pos1,pos2,pos3,pos4},{t,0,2},PlotRange->{0.9,1.1},PlotLegends->{"viscous","Coulomb","Stribeck","static"}]

OpenersExampleSubsection24(Spacer1)Electrical Systems  (6)

Simulate the response of an RLC circuit to a step in the voltage (Subscript[v, 1]) at time t==0.01:

Use component laws together with Kirchhoff's laws for connections:

components={r iR[t]==vR[t],l iL'[t]==vL[t],iC[t]==c vC'[t]};
connections={vR[t]+vL[t]+vC[t]-v1[t]==0,iL[t]==iR[t],iR[t]==iC[t]};

Simulate a step response:

ic={iR[0]==0,iL[0]==0,iC[0]==0,vR[0]==0,vL[0]==0,vC[0]==0};
params={r->5,l->10^-2,c->10^-4};

v1[t_]:=UnitStep[t-.01];

sol=NDSolve[{components,connections,ic}/.params,vC,{t,0,.03}];

Plot[Evaluate[{v1[t],vC[t]}/.sol],{t,0,.03},Filling->{1->0},Ticks->{Range[0.01,0.03,0.01],Automatic},PlotRange->{0,1.5}]

Simulate the response of an RLC circuit to a step in the voltage (Subscript[v, 1]) at time t==0.01:

components2={r1 iR1[t]==vR1[t],l iL'[t]==vL[t],iC[t]==c vC'[t],r2 iR2[t]==vR2[t]};
connections2={vL[t]+vR2[t]==v1[t],vR1[t]+vC[t]==vR2[t],iL[t]==iR1[t]+iR2[t],iR1[t]==iC[t]};

Simulate a step response:

ic2={iR1[0]==0,iR2[0]==0,iL[0]==0,iC[0]==0,vR1[0]==0,vR2[0]==0,vL[0]==0,vC[0]==0};
params2={r1->5,r2->10,l->10^-2,c->10^-4};

v1[t_]:=UnitStep[t-.01];

sol2=NDSolve[{components2,connections2,ic2}/.params2,vC,{t,0,.03}];

Show the step response:

Plot[Evaluate[{v1[t],vC[t]}/.sol2],{t,0,.03},Filling->{1->0},Ticks->{Range[0.01,0.03,0.01],Automatic},PlotRange->{0,1.5}]

Simulate the behavior of a parallel RLC circuit:

components={r iR[t]==vR[t],l iL'[t]==vL[t],iC[t]==c vC'[t]};
connections={iC[t]+iR[t]+iL[t]==i1[t],v[t]==vR[t],vR[t]==vL[t],vL[t]==vC[t]};

Show the response under a constant input current:

i1[t_]:=1;

ic={v[0]==vR[0]==vR[0]==vC[0]==iR[0]==iL[0]==iC[0]==0};
params={r->10,c->10^-3,l->0.4};
sol=NDSolve[{components,connections,iL[0]==0,v[0]==0}/.params,{iR,iL,iC,v},{t,0,0.2},AccuracyGoal->7];

Show the currents in the R, L, C components and the resulting voltage:

{Plot[Evaluate[{iR[t],iL[t],iC[t]}/.sol],{t,0,0.2}],Plot[v[t]/.sol,{t,0,0.2}]}

Model a transistor-amplifier circuit:

The input voltage ve varies sinusoidally:

ve[t_]:=(4/10) Sin[200 Pi t];
params={vb->6,r0->1000,r1->9000,r2->9000,r3->9000,r4->9000,r5->9000,α->99/100,β->10^-6,c1->10^-6,c2->2 10^-6,c3->3 10^-6};

The transistor dispatches the voltage in a nonlinear way, depending on (Subscript[v, 2])-(Subscript[v, 3]):

v23[t_]=β (Exp[(v2[t]-v3[t])/.026]-1);

Use Ohm's law and Kirchhoff's law to determine the governing equations for each node:

node1=c1 (v2'[t]-v1'[t])==v1[t]/r0-ve[t]/r0;
node2=c1(v1'[t]-v2'[t])==(1-α)v23[t]+v2 [t]/r1+v2[t]/r2-vb/r2;
node3=c2 v3'[t]==v23[t]-v3[t]/r3;
node4=c3(v4'[t]-v5'[t])==vb/r4-v4[t]/r4-α v23[t];
node5=c3(v4'[t]-v5'[t])==v5[t]/r5;
ics={v1[0]==0,v2[0]==vb/2,v3[0]==vb/2,v4[0]==vb,v5[0]==0};

Simulate the singular system of equations:

sol=NDSolve[{node1,node2,node3,node4,node5,ics}/.params,{v1,v2,v3,v4,v5},{t,0,0.2},Method->{"EquationSimplification"->"Residual"}] ;

The transistor amplifies voltage (Subscript[v, 4]) relative to (Subscript[v, 3]):

Plot[Evaluate[{v3[t],v4[t]}/.sol],{t,0,0.2},PlotRange->All]

Model a DC-to-DC boost converter from input voltage level vi to desired output voltage level vd using a pulse-width modulated feedback control q[t]:

(RowDefault)

Use Kirchhoff's laws to get a model for the circuit above:

system={vo'[t] ==q[t]i[t]/c -vo[t]/(r c),i'[t]==-q[t]vo[t]/l +vi/l,vo[0]==0,i[0]==0};

The control signal q[t] will switch the transistor on for a fraction vi/vd of each period τ:

control={q[0]==1,WhenEvent[Mod[t,τ]==(vi/vd)τ,q[t]->0],WhenEvent[Mod[t,τ]==0,q[t]->1]};

Boost from a lower voltage vi=24 to a higher voltage vd=36:

pars={vi->24,vd->36,r->22,l->2 10^-1,c->1 10^-4,τ->5 10^-5};

sol=NDSolve[{system,control}/.pars,{vo,i,q},{t,0,.2}, DiscreteVariables->q];

{Plot[{vo[t]/.sol,vd/.pars},{t,0,.2}, PlotRange->All],Plot[q[t]/.sol,{t,0,.0004}]}

Model a DC-to-DC buck-boost converter from input voltage level vi to desired output voltage level vo using a pulse-width modulated feedback control q[t]:

(RowDefault)

Use Kirchhoff's laws to get a model for the circuit above:

system={i'[t]==q[t]vi/l-(1-q[t])vo[t]/l ,vo'[t] ==i[t]/c -vo[t]/(r c),vo[0]==0,i[0]==0};

The control signal q[t] will switch the transistor on for a fraction vd/(vi+vd) of each period:

control={q[0]==1,WhenEvent[Mod[t,τ]==(vd/(vi+vd))τ,q[t]->0],WhenEvent[Mod[t,τ]==0,q[t]->1]};

Boost from a lower voltage vi=24 to a higher voltage vd=36:

pars={vi->24,vd->36,r->22,l->2 10^-2,c->4 10^-6,τ->5 10^-5};

sol=NDSolve[{system,control}/.pars,{vo,i,q},{t,0,.2}, DiscreteVariables->q];

{Plot[{vo[t]/.sol,vd/.pars},{t,0,.02}, PlotRange->All],Plot[q[t]/.sol,{t,0,.0004}]}

Buck from a higher voltage vi=24 to a lower voltage vd=16:

pars={vi->24,vd->16,r->22,l->2 10^-2,c->4 10^-6,τ->5 10^-5};

sol=NDSolve[{system,control}/.pars,{vo,i,q},{t,0,.2}, DiscreteVariables->q];

{Plot[{vo[t]/.sol,vd/.pars},{t,0,.02}, PlotRange->All],Plot[q[t]/.sol,{t,0,.0004}]}

OpenersExampleSubsection25(Spacer1)Hydraulic Systems  (3)

Model the change in height of water in two cylindrical tanks as water flows from one tank to another through a pipe:

Use pressure relations p==ρ g h  and mass conservation:

pressureEqns={p1[t]== ρ g h1[t],p2[t]==ρ g h2[t]};
massConservation={a1 h1'[t]==-flowRate,a2 h2'[t]==flowRate};

Model the flow across the pipe with the Hagen–Poiseuille relation:

flowRate=(p1[t]-p2[t])(π pipeDia^4)/(128 μ pipeLen);

Simulate the system:

params={pipeLen->0.1,pipeDia->0.2,ρ->0.2,μ->2 10^-3,g->9.81};
{a1,a2}={1,1};

sol=NDSolve[{pressureEqns,massConservation,h1[0]==1,h2[0]==0}/.params,{h1,h2,p1,p2},{t,0,15}] ;

Plot[Evaluate[{h1[t],h2[t]}/.sol],{t,0,15}, PlotRange -> All]

Model the second tank as a leaky tank:

massConservation={a1 h1'[t]==-flowRate,a2 h2'[t]==flowRate-a2 h2[t]/10};

sol=NDSolve[{pressureEqns,massConservation,h1[0]==1,h2[0]==0}/.params,{h1,h2,p1,p2},{t,0,80}] ;

Due to the leak in the second tank, both tanks will eventually drain out:

Plot[Evaluate[{h1[t],h2[t]}/.sol],{t,0,80}, PlotRange -> All,AxesOrigin -> {0,0}]

Model the change in height of water in two hemispherical tanks as water flows from one tank to another through a pipe:

pressureEqns={p1[t]== ρ g h1[t],p2[t]==ρ g h2[t]};
massConservation={a1 h1'[t]==-flowRate,a2 h2'[t]==flowRate};

flowRate=(p1[t]-p2[t])(π pipeDia^4)/(128 μ pipeLen);
params={pipeLen->0.1,pipeDia->0.2,ρ->0.2,μ->2 10^-3,g->9.81};

r1=1;r2=1;
hemisphereArea[r_,h_]:=Pi (2 r h-h^2);
a1=hemisphereArea[r1,h1[t]];
a2=hemisphereArea[r2,h2[t]];

sol2=NDSolve[{pressureEqns,massConservation,h1[0]==1,h2[0]==10^-5}/.params,{h1,h2,p1,p2},{t,0,15}] ;

Plot[Evaluate[{h1[t],h2[t]}/.sol2],{t,0,15},PlotRange->All]

Model the change in height of water in three tanks such that one tank feeds water to the other tanks:

pressureEqns={p1[t]==ρ g h1[t],p2[t]==ρ g h2[t],p3[t]==ρ g h3[t]};
massConservation={a1 h1'[t]==-flowRate,a2 h2'[t]==flowRate1,a3 h3'[t]==flowRate2};

Model the flow across the pipe with the Hagen–Poiseuille relation:

flowRate=(p1[t]-pB[t])(π pipeDia^4)/(128μ pipeLenB);
flowRate1=(pB[t]-p2[t])(π pipeDia^4)/(128μ pipeLen1);
flowRate2=(pB[t]-p3[t])(π pipeDia^4)/(128μ pipeLen2);

The flow rate from the first pipe will equal the sum of the flow rates in other two pipes:

flowRateConstraint=flowRate==flowRate1+flowRate2;

Simulate the system:

params={pipeLenB->0.1,pipeLen1->0.1,pipeLen2->0.2,pipeDia->.2,ρ->.2,μ->2 10^-3,g->9.81};
{a1,a2,a3}={1,1,1};

sol=NDSolve[{pressureEqns,massConservation,flowRateConstraint,h1[0]==1,h2[0]==0,h3[0]==0}/.params,{h1,h2,h3},{t,0,15}];

Plot[Evaluate[{h1[t],h2[t],h3[t]}/.sol],{t,0,15},PlotRange->All]

Model the third tank as a leaky tank:

massConservation={a1 h1'[t]==-flowRate,a2 h2'[t]==flowRate1,a3 h3'[t]==flowRate2-a3 h3[t]/10};

sol=NDSolve[{pressureEqns,massConservation,flowRateConstraint,h1[0]==1,h2[0]==0,h3[0]==0}/.params,{h1,h2,h3},{t,0,90}] ;

The first two tanks reach equilibrium and then drain at the same rate:

Plot[Evaluate[{h1[t],h2[t],h3[t]}/.sol],{t,0,90},PlotRange->All]

OpenersExampleSubsection26(Spacer1)Chemical Systems  (2)

Model the kinetics of an autocatalytic reaction:

The rate equations are given as:

{r1,r2,r3}={k1 a[t],k2 b[t]^2,k3 b[t] c[t]};
eqns={a'[t]==-r1+r3,b'[t]==r1-r2-r3};

The concentration of the species a, b, c should always be a constant:

eqEqn={a[t]+b[t]+c[t]==1};
icEqn={a[0]==1,b[0]==0,c[0]==0};

Solve and visualize the evolution of the three species:

params={k1->0.04,k2->3 10^7,k3->10^4};
sol=NDSolve[{eqns,eqEqn,icEqn}/.params,{a,b,c},{t,0,10^8}];

LogLinearPlot[Evaluate[#[t]/.sol],{t,10^-6,10^8}, PlotRange->All,PlotLabel->#[t]]& /@ {a,b,c}

Model a chemical process of two species, FLB and ZHU, that are continuously mixed with carbon dioxide:

The inflow of carbon dioxide per volume unit is denoted by:

Fin=klA (pCO2/H - CO2[t]);

The rate equations are given as:

r1=k1 ((FLB[t])^2)((CO2[t])^0.5);
r2=k2 FLBT[t] ZHU[t];
r3=(k2/KK)FLB[t] ZLA[t];
r4=k3 FLB[t]((ZHU[t])^2) CO2[t];
r5=k4 FLBZHU[t]((CO2[t])^0.5);
eqns={FLB'[t]==-2 r1+r2-r3-r4,CO2'[t]==-0.5 r1-r4-0.5r5+Fin,FLBT'[t]==r1-r2+r3,ZHU'[t]==-r2+r3-2 r4,ZLA'[t]==r2-r3+r5};

Equilibrium equation is given as:

eqEqn=Ks FLB[t] ZHU[t]==FLBZHU[t];

Solve the equations and determine the concentration change in FLB, ZHU, CO2, and ZLA:

params={k1->18.7,k2->0.58,k3->0.09,k4->0.42,KK->34.4,klA->3.3,Ks->115.83,pCO2->0.9,H->737};

ic={FLB[0]==0.444,CO2[0]==0.00123,FLBT[0]==0,ZHU[0]==0.007,ZLA[0]==0};

sol=NDSolve[{eqns,eqEqn,ic}/.params,{FLB,ZHU,,CO2,ZLA},{t,0,200}];

Plot[Evaluate[#[t]/.sol],{t,0,200}, PlotRange->All,PlotLabel->#[t]]& /@ {FLB,ZHU,CO2,ZLA}

OpenersExampleSection5(Spacer1)Properties & Relations  (7)

Symbolic versus numerical differential equation solving:

DSolve[{y''[x]+y[x]==Cos[x],y[0]==y'[0]==0},y,x]

NDSolve[{y''[x]+y[x]==Cos[x],y[0]==y'[0]==0},y,{x,0,20}]

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

The defining equation for (RefLink[JacobiSN, paclet:ref/JacobiSN]):

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

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

Plot[JacobiSN[x,1/2],{x,0,1.5}]

Numerically compute values of an integral at different points in an interval:

iv = Table[{x, NIntegrate[Sin[Exp[s]], {s, 0, x}]}, {x,0,3,.25}]

For functions of the independent variable, (RefLink[NDSolve, paclet:ref/NDSolve]) effectively gives an indefinite integral:

sol = NDSolve[{y'[x] == Sin[Exp[x]],y[0] == 0},  y, {x, 0, 3}]

Plot[Evaluate[y[x] /. sol], {x,0,3}, Epilog->Map[Point, iv]]

Finding an event is related to finding a root of a function of the solution:

ev[t_?NumberQ]:=First[y[t]/.NDSolve[{y''[x]+y[x]==0,y[0]==1,y'[0]==0},y,{x,t}]]

Plot[ev[t],{t,0,3}]

FindRoot[ev[t],{t,1}]

Event location finds the root accurately and efficiently:

NDSolve[{y''[x] + y[x] == 0, y[0] == 1, y'[0] == 0, WhenEvent[y[x]==0, "StopIntegration"]},y,{x, ∞}]

This gives y(10) as a function of y'(0) for a differential equation:

(Subscript[y, 10])[yp0_?NumberQ] := First[y[10] /. NDSolve[{y''[x] + Sin[y[x]] == 0, y[0] == 1, y'[0] == yp0},y,{x,10}]]

Plot[(Subscript[y, 10])[(Subscript[y, 0]^′)], {(Subscript[y, 0]^′), -2,2}]

Find a root of (Subscript[y, 10]) ((Subscript[y, 0])')==0:

FindRoot[(Subscript[y, 10])[(Subscript[y, 0]^′)] == 0, {(Subscript[y, 0]^′),-1.5}]

Solve the equivalent boundary value problem:

Plot[Evaluate[{y[x], y'[x]} /. NDSolve[{y''[x] + Sin[y[x]] == 0, y[0] == 1,y[10] == 0}, y,x]],{x,0,10}]

Use (RefLink[NDSolve, paclet:ref/NDSolve]) as a solver for a (RefLink[SystemModel, paclet:ref/SystemModel]):

sim=SystemModelSimulate[,Method->"NDSolve"]

Plot variables from the simulation result:

SystemModelPlot[sim,{"y","z"}]

Use (RefLink[SystemModel, paclet:ref/SystemModel]) to model larger hierarchical models:

sim=SystemModelSimulate[]

Plot the tank levels in the tank system over time:

SystemModelPlot[sim,{"tank1.h","tank2.h","tank3.h"},PlotRange->All]

OpenersExampleSection6(Spacer1)Possible Issues  (11)

OpenersExampleSubsection27(Spacer1)Numerical Error  (3)

The error tends to grow as you go further from the initial conditions:

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

Find the difference between numerical and exact solutions:

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

Error for a nonlinear equation:

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

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

For high-order methods, the default interpolation may have large errors between steps:

sol = NDSolve[{y''[x]==y[x](y[x]^2-3/2),y[0]==0,y'[0]==1},y,{x,0,10},Method->{"ExplicitRungeKutta", "DifferenceOrder"->8}]

Plot[(y''[x] - y[x](y[x]^2-3/2)) /. sol, {x,0,10}, PlotRange->All]

Interpolation with the order corresponding to the method reduces the error between steps:

dsol = NDSolve[{y''[x]==y[x](y[x]^2-3/2),y[0]==0,y'[0]==1},y,{x,0,10},Method->{"ExplicitRungeKutta", "DifferenceOrder"->8}, InterpolationOrder->All]

Plot[(y''[x] - y[x](y[x]^2-3/2)) /. dsol, {x,0,10}, PlotRange->All]

OpenersExampleSubsection28(Spacer1)Differential Algebraic Equations  (3)

(RefLink[NDSolve, paclet:ref/NDSolve]) cannot automatically handle systems of index greater than 1:

s=NDSolve[{x'''[t]+y[t]==Cos[t],x[t]==Sin[t],x[0]==0,x'[0]==1,x''[0]==0},{x,y},{t,0,2 Pi}]

High-index systems can be solved by performing index reduction on the system:

s=NDSolve[{x'''[t]+y[t]==Cos[t],x[t]==Sin[t],x[0]==0,x'[0]==1,x''[0]==0},{x,y},{t,0,2 Pi},Method->{"IndexReduction"->Automatic}]

Plot[Evaluate[{x[t],y[t]}/.s],{t,0,2 Pi}]

Here is a system of differential-algebraic equations:

DAE = {(Subscript[x, 1])'[t]==(Subscript[x, 3])[t],(Subscript[x, 2])[t] (1-(Subscript[x, 2])[t])==0,(Subscript[x, 1])[t] (Subscript[x, 2])[t]+(Subscript[x, 3])[t] (1-(Subscript[x, 2])[t])==t};

Find the solution with (Subscript[x, 1])'(0)==1:

sol = NDSolve[{DAE,(Subscript[x, 1])'[0]== 1},{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}, {t,0,1}]

(RefLink[NDSolve, paclet:ref/NDSolve]) may change the specified initial conditions if it cannot find the solution with (Subscript[x, 1])'(0)==0:

sol = NDSolve[{DAE,(Subscript[x, 1])'[0]== 0},{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}, {t,0,1}]

(Subscript[x, 1])'[0] /. sol

Change the initial starting guess for the iterations to avoid such issues:

sol = NDSolve[{DAE,(Subscript[x, 1])'[0]==0},{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])},{t,0,1}, Method -> {"DAEInitialization" -> {"Collocation", "DefaultStartingValue" -> 0}}]

(RefLink[NDSolve, paclet:ref/NDSolve]) is limited to index 1, but the solution with (Subscript[x, 2])(0)==1 has index 2:

sol = NDSolve[{DAE,(Subscript[x, 1])[0] == 0, (Subscript[x, 2])[0]== 1, (Subscript[x, 3])[0] == 1},{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}, {t,0,1}]

To solve high-index systems, use index reduction to reduce the DAE to index 1:

sol = NDSolve[{DAE,(Subscript[x, 1])[0] == 0, (Subscript[x, 2])[0]== 1, (Subscript[x, 3])[0] == 1},{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}, {t,0,1}, Method -> {"IndexReduction" -> True,"DAEInitialization" -> "Collocation"}]

The default method may not be able to converge to the default tolerances:

s=100;
dae={x'[t]==-(xs[t]-ys[t]),y'[t]==(xs[t]-ys[t]),xs[t]==s x[t],ys[t]==s y[t],x[0]==1,y[0]==0,xs[0]==s,ys[0]==0};

NDSolve[dae,{x,y},{t,0,1}]

With lower (RefLink[AccuracyGoal, paclet:ref/AccuracyGoal]) and (RefLink[PrecisionGoal, paclet:ref/PrecisionGoal]) settings, a solution is found:

NDSolve[dae,{x,y},{t,0,1},AccuracyGoal->6,PrecisionGoal->6]

The "StateSpace" time integration method can solve this with default tolerances:

NDSolve[dae,{x,y},{t,0,1},Method->{"TimeIntegration"->"StateSpace"}]

OpenersExampleSubsection29(Spacer1)Partial Differential Equations  (3)

Define a nonlinear PDE:

L=20;
wpde={(∂_(t,t))u[t,x]==(∂_(x,x))u[t,x]+(1-((u[t,x])^2)) (1+2 u[t,x]),u[0,x]==(E^(-(x^2))),Derivative[1,0][u][0,x]==0,u[t,-L]==u[t,L]};

The spatial discretization is based on the initial value, which varies less than the final value:

s1=NDSolve[wpde,u,{t,0,L},{x,-L,L}]

By increasing the minimal number of spatial grid points, you can accurately compute the final value:

s2=NDSolve[wpde,u,{t,0,20},{x,-L,L},Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MinPoints"->750}}}]

The plot demonstrates the onset of a spatially more complex solution:

DensityPlot[Evaluate[u[-t,x]/.First[s2]],{x,-L,L},{t,0,-L},PlotPoints->50]

Define a heat equation with an initial value that is a step function:

hpde={D[u[t,x],t] ==  D[u[t,x],x,x], u[0,x] == UnitStep[1/2 - x], u[t,0] == 1, u[t,1] == 0};

Discontinuities in the initial value may result in too many spatial grid points:

Plot3D[Evaluate[u[t,x] /. NDSolve[hpde, u, {t,0,.01},{x,0,1}]],{x,0,1}, {t,0,.01}]

Setting the number of spatial grid points smaller results in essentially as good a solution:

Plot3D[Evaluate[u[t,x] /. NDSolve[hpde, u, {t,0,.01},{x,0,1}, Method->{"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "MaxPoints"->100}}}]],{x,0,1}, {t,0,.01}]

Define a Laplace equation with initial values:

lpde={D[u[x,y],x,x] + D[u[x,y],y,y] == 0, u[x,0] == Sin[Pi x], Derivative[0,1][u][x,0] == 0, u[0,y] == u[1,y] == 0};

The solver only works for equations well posed as initial value (Cauchy) problems:

sol = NDSolve[lpde,u,{x,0,1}, {y,0,1}]

The ill-posedness shows up as numerical instability:

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

OpenersExampleSubsection30(Spacer1)Boundary Value Problems  (1)

This finds a trivial solution of a boundary value problem:

Plot[Evaluate[{y[x], y'[x]} /. NDSolve[{y''[x] + Sin[y[x]] == 0, y[0] == 0,y[10] == 0}, y,x]],{x,0,10}]

You can find other solutions by giving starting conditions for the solution search:

Table[Plot[Evaluate[{y[x], y'[x]} /. NDSolve[{y''[x] + Sin[y[x]] == 0, y[0] == 0,y[10] == 0}, y,x, Method->{"BoundaryValues"->{"Shooting", "StartingInitialConditions"->{y[0] == 0, y'[0] ==i}}}]],{x,0,10}],{i,1,2}]

OpenersExampleSubsection31(Spacer1)Definitions for Unknown Functions  (1)

Definitions for an unknown function may affect the evaluation:

y[x_]:=Cos[x]

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

Clearing the definition for the unknown function fixes the issue:

Clear[y]

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

({{DocumentationFooterIconsRelatedFunction.png, ({{See Also}, {(RefLink[DSolve, paclet:ref/DSolve])  ▪  (RefLink[NDSolveValue, paclet:ref/NDSolveValue])  ▪  (RefLink[ParametricNDSolve, paclet:ref/ParametricNDSolve])  ▪  (RefLink[AsymptoticDSolveValue, paclet:ref/AsymptoticDSolveValue])  ▪  (RefLink[NDEigensystem, paclet:ref/NDEigensystem])  ▪  (RefLink[NDEigenvalues, paclet:ref/NDEigenvalues])  ▪  (RefLink[NIntegrate, paclet:ref/NIntegrate])  ▪  (RefLink[NSolve, paclet:ref/NSolve])  ▪  (RefLink[DifferentialRoot, paclet:ref/DifferentialRoot])  ▪  (RefLink[StreamPlot, paclet:ref/StreamPlot])  ▪  (RefLink[SystemModelSimulate, paclet:ref/SystemModelSimulate])  ▪  (RefLink[ItoProcess, paclet:ref/ItoProcess])}})}})

({{DocumentationFooterIconsRelatedTutorial.png, ({{Tutorials}, {▪(RefLinkPlain[Numerical Differential Equations, paclet:tutorial/IntroductionToNumericalDifferentialEquations])}, {▪(RefLinkPlain[Numerical Solution of Differential Equations, paclet:tutorial/NumericalSolutionOfDifferentialEquations])}, {▪(RefLinkPlain[Advanced Numerical Differential Equation Solving in the Wolfram Language, paclet:tutorial/NDSolveOverview])}, {▪(RefLinkPlain[Numerical Mathematics in the Wolfram Language, paclet:tutorial/NumericalMathematicsInTheWolframLanguage])}})}})

({{DocumentationFooterIconsRelatedGuide.png, ({{Related Guides}, {▪(RefLinkPlain[Differential Equations, paclet:guide/DifferentialEquations])}, {▪(RefLinkPlain[Partial Differential Equations, paclet:guide/PDEModelingAndAnalysis])}, {▪(RefLinkPlain[Equation Solving, paclet:guide/EquationSolving])}, {▪(RefLinkPlain[Solvers over Regions, paclet:guide/GeometricSolvers])}, {▪(RefLinkPlain[Geometric Computation, paclet:guide/GeometricComputation])}, {▪(RefLinkPlain[Mesh-Based Geometric Regions, paclet:guide/MeshRegions])}, {▪(RefLinkPlain[Systems Modeling, paclet:guide/SystemsModeling])}, {▪(RefLinkPlain[Physics & Chemistry: Data and Computation, paclet:guide/PhysicsAndChemistryDataAndComputation])}, {▪(RefLinkPlain[Numerical Evaluation & Precision, paclet:guide/NumericalEvaluationAndPrecision])}, {▪(RefLinkPlain[Calculus, paclet:guide/Calculus])}, {▪(RefLinkPlain[Differential Operators, paclet:guide/DifferentialOperators])}, {▪(RefLinkPlain[Scientific Models, paclet:guide/ScientificModels])}, {▪(RefLinkPlain[Life Sciences & Medicine: Data & Computation, paclet:guide/LifeSciencesAndMedicineDataAndComputation])}})}})

({{(Spacer1), ({{(Spacer1)Related Links}, {(Spacer1)▪(RefLinkPlain[(Spacer1)Implementation notes: Numerical and Related Functions, paclet:tutorial/SomeNotesOnInternalImplementation#3826])}, {(WebLink[DocumentationFooterIconsRelatedLinks-NKS.png(Spacer1)NKS|Online, http://www.wolframscience.com/nksonline/index/search.cgi?SearchIndex=NDSolve]) ((WebLink[A New Kind of Science, http://www.wolframscience.com/nksonline/]))}})}})