Numerical Procedures

William M. MacDonald © 1999

Evaluating Integrals

Real integrals

The numerical integration procedures are very sophisticated and accurate even with the default values for parameters that control the accuracy.

N[NIntegrate[Sqrt[x],{x,0,2}],20]

N[2/3 2^(3/2),20]

Let us try an integral that is singular at the origin.

N[NIntegrate[1.0/Sqrt[x],{x,0,2}],20]

This integral is actually integrable, so Mathematica is not lying.

N[2 2^(1/2),20]

Definite integrals with infinite limits.

Integrals that decay rapidly are easily done.

N[NIntegrate[x^2 Exp[-x^2],{x,0,Infinity}],20]

The exact answer is Sqrt[Pi]/4.

N[Sqrt[Pi]/4,20]

The answer is correct to 8 significant figures.

N[NIntegrate[Log[x] Exp[-x],{x,0,Infinity}],20]

The exact answer is -EulerGamma.

-N[EulerGamma,20]

The answer is correct to 9 significant figures.

A commonly encountered definite integral with an infinite upper limit is the following:

N[NIntegrate[Sin[x]/x,{x,0,Infinity}],20]

The integrand decays too slowly for Mathematica to evalutate. But Mathematica can do some integrals like this exactly by using its tables.

Integrate[Sin[x]/x,{x,0,Infinity}]

Complex integrals

N[NIntegrate[y Exp[I y],{y,0,Pi/4}],20]

The exact integral is (1 - I y) Exp[I y], as Mathematica can tell us.

Integrate[y Exp[I y],y]

Let us check the numerical accuracy by using the exact expression.

N[(Exp[I y](1- I y)/.y->Pi/4) - 
(Exp[I y](1- I y)/.y->0),20]

This looks like a cheat! Let's check with an integrand that has no indefinite integral.

f[y_]=Integrate[Exp[I y-y^2],y]

?Erfi

This is a function that has been studied, but there are no closed forms for it.

N[NIntegrate[Exp[I y - y^2],{y,0,Infinity}],20]

N[f[Infinity]-f[0],20]

This looks like the numerical integration is done approximately rather than using the indefinite integral in terms of the imaginary error function.

Numerical Solution of ODE's

Finding extrema

Finding a minimum

The function FindMinimum is similar to FindRoot.

FindMinimum[4(x^-12 - x^-6),{x,1}]

Plot[4(x^-12 - x^-6),{x, 0.1, 3.0},
PlotRange->{Automatic,{-2.0,2.0}}];

Finding a maximum

Of course, we can use FindMinimum to find maxima also, by looking for the minimum in -f[x].

Clear[x];
FindMinimum[-x Exp[-x],{x,.5}]

Check this. Note that FindMinimum has two components, the value of the function at the minimum and the substitution rule for the value.

D[x Exp[-x],x]/.%[[2]]

Look at the plot.

Plot[x Exp[-x],{x,0,2}]

One interesting feature is that we can give Mathematica fits if we give it a starting value that happens to hit the minimum.

Clear[x];
FindMinimum[-x Exp[-x],{x,1}]

The explanation is that Mathematica first tries to use Newton's method and it fails if the gradient is zero. Obviously Mathematica found a solution to its problem.

Finding multiple minima

Muliple minima  can be found by using the Mathematica function Map with the anonymous function
             FindMinimum[ f[x] , {x, #}]&
as its first argument and a list of estimated roots for the second argument.

Map[ FindMinimum[ Sin[x],{x,#}]&, {5, 8, 11}]

Fitting and interpolating data

Linear least-squares fitting

Here is some data to fit.

dataset=Table[{x,x(1-x)(1.0+Random[Real,{-0.2,0.2}])},
 {x,0,2,0.1}]

g1=ListPlot[dataset,PlotStyle->PointSize[0.02]]

Now fit this data using Fit.

lsqfit=Fit[dataset,{1,x,x^2},x]

g2=Plot[lsqfit,{x,0,2}]

Show[g1,g2]

Interpolation of data

InterpolatingPolynomial

This object will a find a polynomial that runs through all of the points. The data we generated will show the folly of using polynomial fits of high degree.

Clear[x,polyfit];
polyfit[x_]=InterpolatingPolynomial[dataset,x];

g3=Plot[polyfit[x],{x,0,2}]

Show[g1,g3]

Interpolation

This object simply finds an expression that will interpolate between values of a data set.  This is exactly what we want when the data is accurate.  For data with errors this object will not do any smoothing.

intfit=Interpolation[dataset]

g4=Plot[intfit[x],{x,0,2}]

Show[g1,g2,g4]

Here we see that Interpolation really does interpolate between data points. It does not provide a  polynomial "fit".

Nonlinear least-squares fitting

Using FindMinumum

resonance=Table[{x,1.0/(1+5(1-x)^2)(1+Random[Real,{-0.05,0.05}])},
     {x,0,2,0.1}]

gresonancedata=
ListPlot[resonance,PlotStyle->PointSize[0.02],
PlotRange->{Automatic,{0,2}}]

Clear[fitfct];
fitfct[x_,a_,b_,c_,d_]:=a+b/(1 + c(d-x)^2)

fitdeviation[a_,b_,c_,d_]=Apply[Plus,
          Map[(#[[2]]-fitfct[#[[1]],a,b,c,d])^2&,resonance]]

fitparams=FindMinimum[fitdeviation[a,b,c,d],{a,0.},{b,1.},
{c,1.},{d,1.}]

gresonancefit=Plot[
 Evaluate[fitfct[x,a,b,c,d]/.fitparams[[2]]],{x,0,2.}]

Show[gresonancedata,gresonancefit]

Using NonlinearFit

The Package Statistics`NonlinearFit` uses the Levenberg-Marquardt algorithm as the default for performing non-linear fits.

Needs["Statistics`NonlinearFit`"]

nlfit=NonlinearFit[resonance,fitfct[x,a,b,c,d], x,
    { {a , 0.1}, {b , 1.0}, {c, 1.0}, {d, 1.0}}]

nlfitresonance=Plot[
 Evaluate[fitfct[x,a,b,c,d]/.nlfit],{x,0,2.}]

Show[gresonancedata,nlfitresonance]