Numerical Integration Methods.

Kenneth Levasseur
Mathematical Sciences
UMass Lowell
Kenneth_Levasseur@uml.edu

Preliminaries

See http://demonstrations.wolfram.com/ComparingBasicNumericalIntegrationMethods/

a=0;
b=1;
f[x_]:=1/(x+1)^2;
c=(b+a)/2;
tline=(((f[b]-f[a]))/(b-a))(x-a)+f[a];
fp=Plot[f[x],{x,a,b},PlotRange->{0,1.2},AxesOrigin->{0,0},Filling->None];tp=Plot[tline,{x,a,b},PlotRange->{0,1.2},AxesOrigin->{0,0},Filling->Axis,FillingStyle->Opacity[0.3]];Show[{fp,tp}]

What does this image say about the relationship between the  “trapazoid error” and concavity?

Manipulate[a=0;
b=1;
f[x_]:=1/(x+1)^2;
c=(b+a)/2;
line=m (x-c)+f[c];
fp=Plot[f[x],{x,a,b},PlotRange->{0,1.2},AxesOrigin->{0,0},Filling->None];lp=Plot[line,{x,a,b},PlotRange->{0,1.2},AxesOrigin->{0,0},Filling->Axis,FillingStyle->Opacity[0.3]];
Show[{fp,lp}],{{m,0},-3 f'[c],3f'[c]}]

What does this animation say about the relationship between the “midpoint error” and concavity?

Analysis of Integration Errors

Definition of IntegerationError function

IntegrationError[f_,a_,b_]:=Module[{true,m,d,mid,trap},
true=Integrate[f,{x,a,b}]//N;m=(a+b)/2;d=b-a;mid=d (f/.{x->m})//N;trap=d ((f/.{x->a})+(f/.{x->b}))/2//N;{f,a,b,true,mid,trap,true-mid,true-trap}]

heads={"f(x)","a","b","exact","midpoint","trapazoid","mid error","trap.error"};

Pool of test functions

pool={{1/x,1,2.},{E^x,-1,1.},{Cos[x],0.0,Pi/2},{x^3,0,1},{x^4,0,1},{Log[x],1,5},{2^x,0,2},{ArcSin[x],0,1/2},{Sqrt[x],1,4}}

Results

(results=(IntegrationError@@#)&/@pool)//Prepend[#,heads]&

Extract Error Data

data=results//Transpose//Take[#,-2]&//Transpose//N

ListPlot[data,PlotStyle->Red]

Fit data to a line

bf=Fit[data,{x,1},x]

Plot[bf,{x,-0.5,0.5},Epilog->{PointSize[0.02],RGBColor[1,0,0],Point/@data}]

Ratios of errors

Clear[a,b];
data/.{{a_,b_}->{1,b/a}}

Let’s approximate this ratio with 2.    What does this observation about the ratios of errors say in terms of the following quantities?
	I = The exact integral
T= The trapazoid approximation of the integral
M= The midpoint approximation of the integral

Use the relationship you get to find an approximate value of I.

Integration by fitting a quadratic to a function at the endpoints and the midpoint of an interval.

Fit a quadratic to the points {a,f(a)},{((a+b)/2),f(((a+b)/2))},{b,f(b)}.  There is lots of messy algebra, that is somewhat simpler to do by hand, but here is the result.

Clear[f,a,b]
p=InterpolatingPolynomial[{ {a,f(a)},{((a+b)/2),f(((a+b)/2))},{b,f(b)}},x]//Simplify

What is remarkable is that if you integrate this quadratic you get Simpson's Rule.  Notice that -(1/6) (a-b)= ((b-a)/6)

Integrate[p,{x,a,b}]

Simpson’s Rule, as releted to the trapazoid and midpoint rules

For any positive integer n, 
			(Subscript[S, 2n])= (1/3)(Subscript[T, n])+ (2/3)(Subscript[M, n]).

Error Estimates for Trapazoid, Midpoint and Simpson’s rules

The following error estimates appear in virtually every calculus textbook in more or less the same form.  The main variation is that what we call (Subscript[S, 2n]) is occasionally called (Subscript[S, n]), which makes the  estimate for (Subscript[E, S]) look different.

Suppose (Abs[f''(x)])≤K  for  a≤x≤b  then the errors in approximation of (∫_a^b)f(x)dx with n  uniform subintervals using the Trapazoid and Midpoints rules satisfy the bounds
	(Abs[(Subscript[E, T])])≤(((K(b-a))^3)/(12(n^2)))    and   (Abs[(Subscript[E, M])])≤(((K(b-a))^3)/(24(n^2)))

If (Abs[(f^((4)))(x)])≤K'  then (Abs[(Subscript[E, S])]) ≤((K'(((b-a))^5))/(180(n^4))) for even n.

“Hermite Integration?”

Consider approximating a function on each subinterval by the polynomial that agrees with the function and it’s first derivative on the endpoints and the midpoint.     
	6 conditions => the polynomial will have degree 5

Mathematica has a function that does this.

?InterpolatingPolynomial

y=InterpolatingPolynomial[{{{0},0,1},{{Pi/4},(Sqrt[2])/2,(Sqrt[2])/2},{{Pi/2},1,0}},{x}]//N

Plot[{y,Sin[x]},{x,0,Pi/2}]

Plot[{y-Sin[x]},{x,0,Pi/2}]

HermiteIntegrate[f_,{a_,b_},x_,1]:=Module[{y},y=InterpolatingPolynomial[{{{a},f[a],f'[a]},{{(((a+b))/2)},f[(((a+b))/2)],f'[(((a+b))/2)]},{{b},f[b],f'[b]}},{x}];Integrate[y,{x,a,b}]]

Clear[f,a,b]

HermiteIntegrate[f,{a,b},x,1]//Simplify

HermiteIntegrate[Sin,{0,Pi/2},x,1]//N

g=Sin[#^2]&

{HermiteIntegrate[g,{0.,1.},x,1]//N,NIntegrate[g[x],{x,0.,1.}]}

HermiteIntegrate[f_,{a_,b_},x_,n_]:=Module[{d},d=((b-a)/n);Total[HermiteIntegrate[f,{a+(#-1)d,a+# d},x,1]&/@Range[n]]]/;(IntegerQ[n]&&n>1)

HermiteIntegrate[Sin,{0.,Pi/2},x,5]//N

HermiteIntegrate[g,{0.,1.},x,5]-NIntegrate[g[x],{x,0.,1.}]

HermiteIntegrate[f,{a,b},x,4]//Simplify

HermiteIntegrate[f,{0,1},x,4]//Simplify

Improper Integrals

(∫_0^∞)(E^(-(x^2)))\[DifferentialD]x

f[x_]:=(E^(-(x^2)))

fList=Map[f,Range[0,7,0.005]];

Differences[fList,4]/((0.005^4))//Max

Therefore,  let K'=12

g=D[f[x],{x,4}]

Plot[g,{x,0,7},PlotRange->All]

Recall:

If (Abs[(f^((4)))(x)])≤K'  then (Abs[(Subscript[E, S])]) ≤((K'(((b-a))^5))/(180(n^4))) for even n.

2 (10^6) 12  (7^5)/180.

%^(1/4)

M=Map[(7./110)f[#-(7/220)]&,Range[7./110,7.,7./110]]//Total

T=(7./110)(f[0.0]/2+Sum[f[7.k/110],{k,1,109}]+f[7.0]/2)

S=(1/3)T+(2/3)M

NumberForm[S,20]

NumberForm[S,16]

NIntegrate[f[x],{x,0,Infinity}]//NumberForm[#,16]&

N[(Sqrt[π])/2,16]