Numerical Calculus

Overview of capabilities. New technology.

Andrew Moylan

Wolfram Research, Inc.

Plot3D[Sin[x^2+y],{x,0,8},{y,0,8},PlotPoints->75,Filling->None,Mesh->None,AxesLabel->{x,y},ViewAngle->0.48225209790910006`,ViewPoint->{-1.4812854395186446`,-1.938055384130986`,2.3451513330083746`},ViewVertical->{0.`,0.`,1.`}]

NIntegrate[Sin[x^2+y],{x,0,10^4},{y,0,10^4}]

Numerical calculus functions in Mathematica

Numerical calculus functions:

Numerical definite integrals: NIntegrate

Numerical solution of differential equations: NDSolve

Local and global optimization, root-finding: FindMinimum, NMinimize, FindRoot

Representation of approximate functions: InterpolatingFunction

Functions are designed to solve a particular probem, not to implement a specific algorithm: NDSolve, NIntegrate, FindRoot, FindMinimum, ..., not RungeKutta, PredictorCorrector, Gauss, NewtonCotes, Brent, Secant, Newton, ConjugateGradient, ...

Each function contains many different numerical methods. Methods are selected automatically through symbolic analysis of the input, and adaptively during computation.

Detailed control over methods is available through Method options. AccuracyGoal and PrecisionGoal options specify absolute and relative error tolerance.

Output is consistently usable elsewhere in Mathematica.

NIntegrate

One and multiple dimensional definition integration:

NIntegrate[x^x,{x,1,5}]

NIntegrate[((BesselJ[3,y])/(x+1)),{x,0,5},{y,0,5}]

NIntegrate uses adaptive methods with special transformations for singularities:

NIntegrate[((Log[x])/(Sqrt[1-x])),{x,0,1}]

Plot[((Log[x])/(Sqrt[1-x])),{x,0,1}]

Integrands can be vector or tensor valued:

NIntegrate[(({{x, (Sqrt[x]), (x^(1/3))}, {(Sqrt[x]), (x^(1/3)), (x^(1/4))}, {(x^(1/3)), (x^(1/4)), (x^(1/5))}})),{x,0,5}]//MatrixForm

Specialized quadrature rules might be automatically selected, depending on the form on the integrand. Symbolic preprocessing also helps in other ways.

One important case is integrands containing discontinuous functions:

Plot[Sin[x^2+Floor[x]],{x,0,5},ExclusionsStyle->{Dashed,Red}]

Discontinuities inside the integration region drastically slow the convergence of standard quadrature rules. NIntegrate automatically subdivides integrals containing piecewise functions:

NIntegrate[Sin[x^2+Floor[x]],{x,0,5}]

Sets of points satisfying an equality can be manually excluded:

Plot[{Re[(1/(Sqrt[Sin[x]]))],Im[(1/(Sqrt[Sin[x]]))]},{x,0,10},Filling->Axis]

NIntegrate[(1/(Sqrt[Sin[x]])),{x,0,10},Exclusions->Sin[x]==0]

NIntegrate: Integration over implicit regions

Using Boole, the integration region can be defined implicitly in terms of a predicate:

Clear[x,y,z]
RegionPlot[x^2+y^2≤1∧y≥0∧y≥-x,{x,-1,1},{y,-1,1},ImageSize->Small]

The integral of a Boole gives the area:

NIntegrate[Boole[x^2+y^2≤1∧y≥0∧y≥-x],{x,-∞,∞},{y,-∞,∞}]

Show where NIntegrate samples the integrand:

Show[{,ListPlot[Last[Reap[NIntegrate[Boole[x^2+y^2≤1∧y≥0∧y≥-x],{x,-1,1},{y,-1,1},EvaluationMonitor:>Sow[{x,y}],Method->{"SymbolicPreprocessing","UnitCubeRescaling"->False}]]],PlotRange->{{-1,1},{-1,1}},AspectRatio->1]}]

Integral of a scalar field over the same region:

NIntegrate[(1+Sin[(x+2)(y+2)])Boole[x^2+y^2≤1∧y≥0∧y≥-x],{x,-1,1},{y,-1,1}]

Plot3D[1+Sin[(x+2)(y+2)],{x,-1,1},{y,-1,1},RegionFunction->({x,y}|->x^2+y^2≤1∧y≥0∧y≥-x),Filling->0]

A 3D example, half of an ellipsoidal shell:

RegionPlot3D[1≤x^2+y^2+1/2z^2≤3∧y≥0,{x,-3,3},{y,-3,3},{z,-3,3},Mesh->None,PlotPoints->100]

Volume:

NIntegrate[Boole[1≤x^2+y^2+1/2z^2≤3∧y≥0],{x,-3,3},{y,-3,3},{z,-3,3}]

Homework

"The base of a solid is the region bounded by the curves y=2sin[x] and y=tan[x] in the first quadrant between A(0,0) and B(pi/3, sqrt[3). The solid is formed such that every cross-sectional slice perpendicular to the x-axis is a semicircle having a diameter with one end-point on the curve y=2sin[x] and the other on y=tan[x]. Find the volume of this solid."

Plot[{2Sin[x],Tan[x]},{x,0,π/3}]

Δy=y-((2Sin[x]+Tan[x])/2);
r=((2Sin[x]-Tan[x])/2);
NIntegrate[Boole[Δy^2+z^2<r^2],{x,0,(π/3)},{y,Tan[x],2Sin[x]},{z,0,∞}]

The region looks like this:

RegionPlot3D[Tan[x]<y<2Sin[x]∧(y-((2Sin[x]+Tan[x])/2))^2+z^2<(((((2Sin[x]-Tan[x])/2)))^2),{x,0,(π/3)},{y,0,2},{z,0,0.3},AxesLabel->{x,y,z},PlotPoints->50,BoxRatios->Automatic]

NIntegrate: Parametric surfaces

A trefoil knot as a parametric curve:

{x,y,z}={(2+Cos[(3/2)t])Cos[t],(2+Cos[(3/2)t])Sin[t],Sin[(3/2)t]};

ParametricPlot3D[{x,y,z},{t,0,4π}]

Arc length:

NIntegrate[(Sqrt[D[x,t]^2+D[y,t]^2+D[z,t]^2]),{t,0,4π}]

General parametric two-dimensional surface embedded in three dimensions:

{x,y,z}={(2+Cos[v])Cos[u],(2+Cos[v])Sin[u],Sin[v](5/4-Cos[u])};

ParametricPlot3D[{x,y,z},{u,0,2π},{v,0,2π}]

The surface area is:

∯ (Norm[((∂{x,y,z})/(∂u))×((∂{x,y,z})/(∂v))]) dudv

This translates straight into Mathematica:

NIntegrate[Norm[Cross[D[{x,y,z},u],D[{x,y,z},v]]],{u,0,2π},{v,0,2π}]

Or like this:

NIntegrate[(Norm[((∂{x,y,z})/(∂u))\[Cross]((∂{x,y,z})/(∂v))]),{u,0,2π},{v,0,2π}]

Stokes theorem lets us work out more properties of the region enclosed by a parametric surface. Reminder: In 3D, it is

∫∫∫∇.Fdxdydz==∯F.((∂{x,y,z})/(∂u))×((∂{x,y,z})/(∂v))dudv

F is any vector field we choose.

To find the volume inside the region, we want ∇.F==1. F=={x,y,z}/3 is suitable:

m=NIntegrate[(({x,y,z})/3).Cross[D[{x,y,z},u],D[{x,y,z},v]],{u,0,2π},{v,0,2π}]

To work out the center of mass we want ∇.F == {x,y,z} (so F should be a matrix). One option is:

cm=(1/m)NIntegrate[DiagonalMatrix[(({x,y,z}^2)/2)].Cross[D[{x,y,z},u],D[{x,y,z},v]],{u,0,2π},{v,0,2π},AccuracyGoal->4]

Show the center of mass:

Show[ParametricPlot3D[{x,y,z},{u,0,2π},{v,0,2π},PlotStyle->Opacity[1/2],Mesh->None],Graphics3D@Line[Transpose[{cm+10IdentityMatrix[3]/3,cm-10IdentityMatrix[3]},{2,3,1}]]]

Clear[x,y,z]

NDSolve: ODEs

See also: Mark Sofroniou's talk specifically on NDSolve: tomorrow 11:30am.

An ordinary differential equation, first order, nonlinear:

sol=y/.First[NDSolve[{y'[x]==x^2+y[x]^2-1,y[0]==0.5},y,{x,-2,2}]]

NDSolve gives the approximate solution function as an InterpolatingFunction, a piecewise polynomial. It works like any other function in Mathematica:

Plot[sol[x],{x,-2,2}]

We can use VectorPlot to composite this with a direction field:

Show[VectorPlot[Normalize[{1,x^2+y^2-1}],{x,-2,2},{y,-2,2}],ParametricPlot[{x,sol[x]},{x,-2,2},PlotStyle->{Thick,Red}]]

StreamPlot, new in Mathematica 7, uses NDSolve behind the scenes:

StreamPlot[{1,x^2+y^2-1},{x,-2,2},{y,-2,2},StreamPoints->{{{{1,0},Red},Automatic}}]

NDSolve: PDEs

A fun example of a partial differential equation is the Schroedinger equation in one dimension:

Clear[x,t]
equations={I D[ψ[x,t],t]==-(1/2)D[ψ[x,t],{x,2}]+250Sin[x/2]^10ψ[x,t],ψ[x,0]==Exp[I 20 x]Exp[-10x^2],ψ[-π,t]==ψ[π,t]};

The differential equation models a particle rolling back and forth in a periodic well shaped like this:

Plot[250Sin[x/2]^10,{x,-π,π},ImageSize->Small]

To solve this partial differential equation NDSolve first discretizes in the spatial direction, converting the problem into a set of coupled ordinary differential equations in the time dimension.

sol=ψ/.First[NDSolve[equations,ψ,{x,-π,π},{t,0,1/2}]]

The solution as a function of position and time:

Plot3D[Evaluate[Abs[sol[x,t]^2]],{x,-π,π},{t,0,1/2}, PlotRange->All,PlotPoints->50,AxesLabel->{x,t}]

An animation is a nice way to view the time evolution of the system:

Animate[Plot[{Sin[(x/2)]^10+Abs[sol[x,t]]^2,Sin[(x/2)]^10},{x,-π,π},PlotRange->{0,3/2},Filling->1->{2}],{t,0,1/2},DefaultDuration->10]

NDSolve: Other types of problems

NDSolve can handle special types of differential equations including differental-algebraic equations and delay differential equations.

Methods for delay differential equations were introduced in Mathematica 7.

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

Plot[Evaluate[{x[t],x'[t]} /. First[sol]],{t,0,5}]

Take a look at the NDSolve advanced documentation for nice examples of DDEs.

There are also algorithms suited to particular kinds of numerically difficult problems. An example for stiff equations:

sol={u[t],v[t]}/.First[NDSolve[{u'[t]==v[t],((3 v'[t])/1000)==-u[t]+(1-((u[t])^2)) v[t],u[0]==2,v[0]==0},{u,v},{t,0,5/2},Method->"StiffnessSwitching"]]

Plot[sol,{t,0,5/2},PlotRange->15]

Mark will be talking about some new methods for stiff problems in his NDSolve talk tomorrow.

NDSolve: Indefinite integration

An indefinite integral is just a very simple initial value problem:

NDSolve[{y'[x]==Exp[Cos[x]],y[0]==0},y,{x,0,10}]

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

This operation can be wrapped up into a function:

NIndefiniteIntegral[expr_,{x_,a_,b_}]:=Module[{y},y/.First[NDSolve[{y'[x]==expr,y[a]==0},y,{x,a,b}]]]

f=NIndefiniteIntegral[Exp[Sin[x]],{x,0,10}]

f is like any other function:

Plot[f[x],{x,0,10}]

FindRoot[f[x]==7,{x,0,10}]

Approximate functions

The InterpolatingFunction objects produced by NDSolve represent approximate functions usually represented internally as piece polynomials. You can also generate them directly with FunctionInterpolation.

Clear[x,y,z]
nf=FunctionInterpolation[Exp[BesselJ[0,x]],{x,0,10}]

Plot[nf[x],{x,0,10}]

You use them like any other function in Mathematica:

{D[nf[x],x],nf[x],Integrate[nf[x],x]}

Plot[%,{x,0,10}]

Integrate[nf[x],{x,0,10}]

NIntegrate[Exp[BesselJ[0,x]],{x,0,10}]

Optimization and root-finding

We can use FindMinimum to demonstrate the integration of several numerical calculus functions.

First define a function that gives the solution to a differential equation that depends on a parameter a:

Clear[x,y,sol]
sol[a_?NumberQ]:=y/.First[NDSolve[{y''[x]==-y[x](1+Cos[x+y[x]]),y[0]==1,y'[0]==a},y,{x,0,10}]]

Consider the solutions for parameter values 0<a<1:

sol[0.5]

Animate[Plot[Evaluate[sol[a][x]],{x,0,10},PlotRange->5],{a,0,1}]

We can numerically explore the properties of these solutions.

With NIntegrate we can compute the arc length of any of these solutions:

arclength[a_?NumericQ]:=NIntegrate[(Sqrt[1+sol[a]'[x]^2]),{x,0,10}]

Quick plot of the arc length as a function of a:

plot=ListLinePlot[Table[{a,arclength[a]},{a,0,1,0.05}]]

The minimum near x=0.8:

min=FindMinimum[arclength[a],{a,0.8,0,1}]

FindMinimum also does constrained local minimization. Multidimensional example:

{z,xy}=FindMinimum[{Sin[x]Sin[2y+1],x^2+y^2<3},{{x,2},{y,2}}]

p=Show[Plot3D[Sin[x]Sin[2y+1],{x,-2,2},{y,-2,2},RegionFunction->({x,y}|->x^2+y^2<3),MeshFunctions->{#3&},AxesLabel->{x,y}],Graphics3D[{Red,PointSize[0.05],Point[{x,y,z}/.xy]}]]

NMinimize does global optimization.

NMinimize[{Sin[x]Sin[2y+1],x^2+y^2<3},{x,y}]

FindRoot performs a local root search. We could get a solution to the differential equation above whose arc length is exactly 15:

FindRoot[arclength[a]==15,{a,0.5}]

Plot[Evaluate[sol[a][x]/.%],{x,0,10}]

FindRoot can be useful for numerically inverting functions that have no standard inverse:

Solve[x^x+x==y,x]

Plot[x^x+x,{x,1,2}]

inverse[y_?NumericQ]:=x/.FindRoot[x^x+x==y,{x,1}]

Plot[inverse[y],{y,2,6}]

New technology: General oscillatory integration

Integrals of highly oscillatory functions are very expensive to compute with traditional quadrature rules.

NIntegrate[(1/(1+x^3))Sin[x],{x,0,∞},Method->"GaussKronrodRule"]

Since Mathematica 6, NIntegrate has had algorithms for various specialized types of oscillatory integrals. Perhaps the most commonly used are those for Fourier integrals. Such a rule is automatically selected here:

NIntegrate[(1/(1+x^3))Sin[x],{x,0,∞}]

New methods we are now developing will allow efficient automatic integration of a great variety of oscillatory functions for the very first time.

First of Trefethen's "SIAM 100 digit challenge" problems:

NIntegrate[(1/x)Cos[((Log[x])/x)],{x,0,1}]

Plot[(1/x)Cos[((Log[x])/x)],{x,0,1}]

Multivariate highly oscillatory integration:

NIntegrate[Sin[Exp[x]+y^2],{x,0,∞},{y,0,∞}]

Plot3D[Sin[Exp[x]+y^2],{x,0,4},{y,0,4},PlotPoints->50]

Oscillatory integration: Levin's method

NIntegrate's new oscillatory integration methods are based on "Levin's method". They are a good example of the kind of methods that are only properly realizable in a system in which numerics and symbolics are unified.

Most oscillatory special functions satisfy a linear differential equation with non-highly-oscillatory coefficients. Often the coeffients are polynomial:

DifferentialRootReduce[Sin[x],x]

DifferentialRootReduce[StruveH[2,x],x]

Levin's method applies to any any sums or products of such functions, as well as their products and compositions with non-oscillatory functions:

NIntegrate[Sin[x] (Sin[x^2] (Sin[x^3]+(1/x))+(1/x)),{x,0,∞}]

NIntegrate[BesselY[2/3,x]Sin[x+(Sqrt[x])]^4,{x,0,10^4}]

These are diffcult functions:

Table[Plot[f,{x,0,20},PlotLabel->f,ImageSize->Medium],{f,{BesselY[2/3,x]Sin[x+(Sqrt[x])]^4,Sin[x] ((1/x)+Sin[x^2] ((1/x)+Sin[x^3]))}}]

To be practically useful, it is essential to automatically determine the linear differential relation satisfied by the integrand. Mathematica is of course uniquely well equipped to do this.

For example, the differential relation found by NIntegrate for sin(x) ((sin((x^3))+(1/x)) sin((x^2))+(1/x)) is:

D[(({{Sin[x] Sin[(x^2)] Sin[(x^3)]}, {Cos[(x^3)] Sin[x] Sin[(x^2)]}, {Cos[(x^2)] Sin[x] Sin[(x^3)]}, {Cos[(x^2)] Cos[(x^3)] Sin[x]}, {Sin[x] Sin[(x^2)]}, {Cos[(x^2)] Sin[x]}, {Cos[x] Sin[(x^2)] Sin[(x^3)]}, {Cos[x] Cos[(x^3)] Sin[(x^2)]}, {Cos[x] Cos[(x^2)] Sin[(x^3)]}, {Cos[x] Cos[(x^2)] Cos[(x^3)]}, {Cos[x] Sin[(x^2)]}, {Cos[x] Cos[(x^2)]}, {Sin[x]}, {Cos[x]}})),x]==(({{0, 3 (x^2), 2 x, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0}, {-3 (x^2), 0, 0, 2 x, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0}, {-2 x, 0, 0, 3 (x^2), 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, {0, -2 x, -3 (x^2), 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 2 x, 0, 0, 0, 0, 1, 0, 0, 0}, {0, 0, 0, 0, -2 x, 0, 0, 0, 0, 0, 0, 1, 0, 0}, {-1, 0, 0, 0, 0, 0, 0, 3 (x^2), 2 x, 0, 0, 0, 0, 0}, {0, -1, 0, 0, 0, 0, -3 (x^2), 0, 0, 2 x, 0, 0, 0, 0}, {0, 0, -1, 0, 0, 0, -2 x, 0, 0, 3 (x^2), 0, 0, 0, 0}, {0, 0, 0, -1, 0, 0, 0, -2 x, -3 (x^2), 0, 0, 0, 0, 0}, {0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 2 x, 0, 0}, {0, 0, 0, 0, 0, -1, 0, 0, 0, 0, -2 x, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0}})).(({{Sin[x] Sin[(x^2)] Sin[(x^3)]}, {Cos[(x^3)] Sin[x] Sin[(x^2)]}, {Cos[(x^2)] Sin[x] Sin[(x^3)]}, {Cos[(x^2)] Cos[(x^3)] Sin[x]}, {Sin[x] Sin[(x^2)]}, {Cos[(x^2)] Sin[x]}, {Cos[x] Sin[(x^2)] Sin[(x^3)]}, {Cos[x] Cos[(x^3)] Sin[(x^2)]}, {Cos[x] Cos[(x^2)] Sin[(x^3)]}, {Cos[x] Cos[(x^2)] Cos[(x^3)]}, {Cos[x] Sin[(x^2)]}, {Cos[x] Cos[(x^2)]}, {Sin[x]}, {Cos[x]}}))

where the integrand is

Sin[x] (Sin[x^2] (Sin[x^3]+(1/x))+(1/x))==(({{1}, {0}, {0}, {0}, {1/x}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {1/x}, {0}})).(({{Sin[x] Sin[(x^2)] Sin[(x^3)]}, {Cos[(x^3)] Sin[x] Sin[(x^2)]}, {Cos[(x^2)] Sin[x] Sin[(x^3)]}, {Cos[(x^2)] Cos[(x^3)] Sin[x]}, {Sin[x] Sin[(x^2)]}, {Cos[(x^2)] Sin[x]}, {Cos[x] Sin[(x^2)] Sin[(x^3)]}, {Cos[x] Cos[(x^3)] Sin[(x^2)]}, {Cos[x] Cos[(x^2)] Sin[(x^3)]}, {Cos[x] Cos[(x^2)] Cos[(x^3)]}, {Cos[x] Sin[(x^2)]}, {Cos[x] Cos[(x^2)]}, {Sin[x]}, {Cos[x]}}))

There are a whole host of integrals involving special functions which can now be done automatically and efficiently:

NIntegrate[((BesselJ[0,x])/(1+x)),{x,0,10^4}]

NIntegrate[Sinh[x],{x,5,10^6I}]

NIntegrate[SinIntegral[x],{x,0,10^4}]

Automatic multidimensional oscillatory integration is also a first:

NIntegrate[Exp[I(x+y)],{x,1,10^4},{y,1,10^4}]

The decision to use oscillatory methods takes into account the region of the complex plane over which a function is being evaluated.

In this example, Power is automatically identified as a highly oscillatory function over the given range:

NIntegrate[2^(x y),{x,I,10^4I},{y,1,10^4}]

In this example, Power is identified as non-oscillatory:

NIntegrate[2^(x y),{x,I,10^4I},{y,I,2I}]

Summary

Integrated design makes it easy to use the output of one function as the input to another.

Problem-based functions that implement a lot of algorithms and switch between them adaptively and through symbolic preprocessing.

Symbolic and numerical methods together go much further than either one by itself.