Math 4401 Numerical Methods

Notebook Author/Course Instructor: Barry McQuarrie
Date: Oct 2016
Lecture topic: Romberg Integration and Adaptive Quadrature

Romberg Integration

Clear[f,R,Rom,b,a,num]

f[x_]:=Sqrt[1-Cos[x]^2]/Sqrt[1+x^3];
Integrate[f[x],{x,0,2}]

Composite Quadrature (trapezoidal rule) to estimate the Integral

trap[f_,a_,b_]:=Module[{h,f0,f1,f2,f3,f4},
h=b-a;
f0=f/.x->a;
f1=f/.x->a+h;
Return[SetPrecision[(1/2) h f0+(1/2) h f1,12]];
]

trapcomposite[f_,a_,b_,n_]:=Module[{h,sum,x1,x2},
h=N[((b-a)/n)];
sum=0.0;
x1=a-h;
x2=a;
Do[
x1=x1+h;
x2=x2+h;
sum=sum+trap[f,x1,x2]
,{i,1,n}];
Return[SetPrecision[sum,12]];
]

TableForm[Table[{n,trapcomposite[f[x],0,2,n]},{n,50,500,50}],TableHeadings->{None,{"Subintervals","Trapezoidal Approximation"}}]

errorTrapezoidal=Table[{n,Abs[trapcomposite[f[x],0,2,n]-trapcomposite[f[x],0,2,n-1]]},{n,10,500,10}];
trapplot=ListLogPlot[errorTrapezoidal,PlotStyle->{Red,PointSize[0.02`]},Frame->True,FrameLabel->{"n","Log Error"},ImageSize->{500,Automatic},BaseStyle->20]

Composite Quadrature (n=4) to estimate the Integral

n4[f_,a_,b_]:=Module[{h,f0,f1,f2,f3,f4},
h=((b-a)/4);
f0=f/.x->a;
f1=f/.x->a+h;
f2=f/.x->a+2h;
f3=f/.x->a+3h;
f4=f/.x->a+4h;
Return[SetPrecision[(14/45) h f0+(64/45) h f1+(8/15) h f2+(64/45) h f3+(14/45) h f4,12]];
]

n4composite[f_,a_,b_,n_]:=Module[{h,sum,x1,x2},
h=N[((b-a)/n)];
sum=0.0;
x1=a-h;
x2=a;
Do[
x1=x1+h;
x2=x2+h;
sum=sum+n4[f,x1,x2]
,{i,1,n}];
Return[SetPrecision[sum,20]];
]

TableForm[Table[{n,n4composite[f[x],0,2,n]},{n,10,50,5}],TableHeadings->{None,{"Subintervals","Trapezoidal Approximation"}}]

result=n4composite[f[x],0,2,50]

errorn4=Table[{n,Abs[n4composite[f[x],0,2,n]-n4composite[f[x],0,2,n-1]]},{n,2,40,2}];
n4plot=ListLogPlot[errorn4,PlotStyle->{Blue,PointSize[0.02`]},Frame->True,FrameLabel->{"n","Log Error"},ImageSize->{500,Automatic},BaseStyle->20]

Romberg Integration

Clear[f,R,Rom,b,a,num]

f[x_]:=Sqrt[1-Cos[x]^2]/Sqrt[1+x^3];
a=0.0;
b=2.0;
num=7;

h=b-a;
R[1,1]=Rom[1,1]=trapcomposite[f[x],a,b,1];

i=2;
Print["Romberg Integration"]
Print["Input:  ","f(x) = ",f[x]," a = ",a," b = ",b," n = ",num]
While[i≤num ,

R[2,1]=trapcomposite[f[x],a,b,(2^(i-1))];

Do[R[2,j]=R[2,j-1]+((R[2,j-1]-R[1,j-1])/(4^(j-1)-1)),{j,2,i}];

h=h/2;

Do[Rom[i,j]=R[2,j],{j,1,i}];
Do[R[1,j]=R[2,j],{j,1,i}];

i=i+1]

data=Table[Rom[i,j],{i,1,num},{j,1,i}];
TableForm[SetPrecision[data,7],TableSpacing->1]

SetPrecision[Rom[num,num],12]

errorrom=Table[{n,Abs[Rom[n,n]-result]},{n,1,num}];
romplot=ListLogPlot[errorrom,PlotStyle->{Purple,PointSize[0.02`]},Frame->True,FrameLabel->{"n","Log Error"},ImageSize->{500,Automatic},BaseStyle->20]

Error Comparison

The n's aren't really the same, but we should still be impressed by how much the Romberg technique speeds up the convergence. But it's not significanty better than the n=4 composite technique.

Show[trapplot,n4plot,romplot,Axes->False,PlotRange->{Automatic,Automatic}]

Show[n4plot,romplot,Axes->False,PlotRange->{Automatic,Automatic}]

Adaptive Trapezoidal Integration

Here is a function where it would be beneficial to put more nodes in regions where the function is changing rapidly, and fewer nodes in regions where the function is not changing rapidly.

f[x_]=Cos[Exp[x]];
Plot[f[x],{x,0,Pi},PlotStyle->{Blue,Thick}]
NIntegrate[f[x],{x,0,Pi}]

Here is our trapezoidal module:

trap[f_,a_,b_]:=Module[{h,f0,f1,f2,f3,f4},
h=b-a;
f0=f/.x->a;
f1=f/.x->a+h;
Return[SetPrecision[(1/2) h f0+(1/2) h f1,12]];
]

Here is a dangerous bit of code, but it should work if your function is not crazy. The code is dangerous because it calls itself and has no checks or breaks in case things go wrong. But, it works for our example, and is simple to understand! Note for each halving the error tolerance goes down by a factor of 2, but the trapezoidal rule error should go down by a factor of (2^3)=8, so after some amount of halvings we should meet the required tolerance.

AdaptiveTrap[f_,a_,b_,tol_]:=Module[{c,Sab,Sac,Sbc},
c=N[(a+b)/2];
Sab=trap[f,a,b];
Sac=trap[f,a,c];
Scb=trap[f,c,b];
If[
Abs[Sab-Sac-Scb]<tol,
Return[Sac+Scb],
Return[AdaptiveTrap[f,a,c,tol/2]+AdaptiveTrap[f,c,b,tol/2]]
];
]

AdaptiveTrap[f[x],0,Pi,0.001]

Let’s modify it to keep track of the c,f(c) values being used. This will be informative.

nodelist={};
functionevallist={};
AdaptiveTrap[f_,a_,b_,tol_]:=Module[{c,Sab,Sac,Sbc},
c=N[(a+b)/2];
nodelist=Append[nodelist,{c,0}];
functionevallist=Append[functionevallist,{c,f/.x->c}];
Sab=trap[f,a,b];
Sac=trap[f,a,c];
Scb=trap[f,c,b];
If[
Abs[Sab-Sac-Scb]<tol,
Return[Sac+Scb],
Return[AdaptiveTrap[f,a,c,tol/2]+AdaptiveTrap[f,c,b,tol/2]]
];
]

AdaptiveTrap[f[x],0,Pi,0.005]

ListPlot[{functionevallist,nodelist},ImageSize->{700,Automatic},PlotStyle->PointSize[0.01],BaseStyle->20]

ListPlot[{functionevallist,nodelist},ImageSize->{700,Automatic},PlotRange->{{2.5,2.8},All},PlotStyle->PointSize[0.015],BaseStyle->20]

Notice where the curve is essentially straight (which trapezoid does very well) there are fewer functional evaluations. Cool.

Adaptive Simpson’s 3/8 Integration

Here is a function where it would be beneficial to put more nodes in regions where the function is changing rapidly, and fewer nodes in regions where the function is not changing rapidly.

f[x_]=Cos[Exp[x]];
Plot[f[x],{x,0,Pi},PlotStyle->{Blue,Thick}]
NIntegrate[f[x],{x,0,Pi}]

Here is our Simpson’s 3/8 module:

simp38[f_,a_,b_]:=Module[{h,f0,f1,f2,f3},
h=((b-a)/3);
f0=f/.x->a;
f1=f/.x->a+h;
f2=f/.x->a+2h;
f3=f/.x->a+3h;
Return[SetPrecision[(3/8) h f0+(9/8) h f1+(9/8) h f2+(3/8) h f3,12]];
]

This should be an improvement over the adaptive trapezoidal, since each halving the error is Simpson’s 3/8 goes down by a factor of (2^5)=32, whereas trapezoidal was only (2^3)=8.

AdaptiveSimp38[f_,a_,b_,tol_]:=Module[{c,Sab,Sac,Sbc},
c=N[(a+b)/2];
Sab=simp38[f,a,b];
Sac=simp38[f,a,c];
Scb=simp38[f,c,b];
If[
Abs[Sab-Sac-Scb]<tol,
Return[Sac+Scb],
Return[AdaptiveSimp38[f,a,c,tol/2]+AdaptiveSimp38[f,c,b,tol/2]]
];
]

And it does seem to be faster.

Timing[AdaptiveSimp38[f[x],0,Pi,(10.0^(-5))]]
Timing[AdaptiveTrap[f[x],0,Pi,(10.0^(-5))]]

Less Dangerous Module Using Simspon’s Rule

This is a straight adaptation of an algorithm for adaprive quadrature I found in the book Numerical Analysis by Burden, Faires. I have a link to a copy of the algorithm on the moodle site.

AdaptiveQuadratureSimpsons[f_,a0_,b0_,tol0_,NL_]:=Module[{i,APP,tol,a,h,FA,FC,FB,S,L,FD,FE,S1,S2,v1,v2,v3,v4,v5,v6,v7,v8},
i=1;
APP=0.0;
tol[i]=10*tol0;
a[i]=N[a0];
h[i]=N[(b0-a0)/2];
FA[i]=N[f/.t->a0];
FC[i]=N[f/.t->(a0+h[i])];
FB[i]=N[f/.t->b0];
S[i]=h[i](FA[i]+4FC[i]+FB[i])/3; (*Approximation over entire interval*)
L[i]=1;

While[i>0,
FD=f/.t->(a[i]+h[i]/2);
FE=f/.t->(a[i]+3h[i]/2);
S1=h[i](FA[i]+4FD+FC[i])/6; (*Approximation over halves of subinterval*)
S2=h[i](FC[i]+4FE+FB[i])/6;
v1=a[i];
v2=FA[i];
v3=FC[i];
v4=FB[i];
v5=h[i];
v6=tol[i];
v7=S[i];
v8=L[i];

i=i-1;(*Delete the level*)

If[Abs[S1+S2-v7]<v6,
APP=APP+(S1+S2),
If[v8≥NL,
Print["Output Level Exceeded"];
Abort[],
i=i+1;(*Data for right half subinterval*)
a[i]=v1+v5;
FA[i]=v3;
FC[i]=FE;
FB[i]=v4;
h[i]=v5/2;
tol[i]=v6/2;
S[i]=S2;
L[i]=v8+1;

i=i+1;(*Data for left half subinterval*)
a[i]=v1;
FA[i]=v2;
FC[i]=FD;
FB[i]=v3;
h[i]=h[i-1];
tol[i]=tol[i-1];
S[i]=S1;
L[i]=L[i-1];
];
];
];
Return[SetPrecision[APP,16]];
]

Let’s check it out!

f[t_]=Cos[3(t^2)]+t;
a=0;
b=2Pi;
Plot[f[t],{t,a,b},PlotStyle->{Blue,Thick}]

Timing[AdaptiveQuadratureSimpsons[f[t],a,b,(10^(-14)),200]]
Timing[SetPrecision[NIntegrate[f[t],{t,a,b}],16]]

Note you should check to see if the routine aborts gracefully:

AdaptiveQuadratureSimpsons[f[t],a,b,(10^(-14)),10]

Be careful with adaptive routines, since they may exit prematurely in some cases. I am going to ask you tio investigate this on the midterm (I think).