5.9 Approximate Integration

Mathematica script by Chris Parrish,
   cparrish@sewanee.edu

Sources and references for some of these problems include
    James Stewart, "Calculus: Concepts and Contexts," Second Edition, Brooks/Cole, 2001
    Deborah Hughes-Hallett, Andrew M. Gleason, et. al., "Calculus," Second Edition, John Wiley & Sons, 1998
    Robert Fraga, ed., "Calculus Problems for a New Century," The Mathematical Association of America, 1993

Mathematica Code for Some
Composite Numerical Integration Formulas (L, R, T, M, S)

Composite Rules -- This section contains Mathematica code only, not worked examples

Using Composite Numerical Integration Formulas (L, R, T, M)

Estimating RowBox[{RowBox[{∫_0, ^, 2.5}], Sin[t^2] x}]

Hughes-Hallett, Gleason, et al, First Edition, Text Example, page 380

Let's approximate

RowBox[{(*,  , RowBox[{RowBox[{∫_0, ^, 2.5}], Sin[t^2] x}],  , *)}]

Clear[f, t, a, b] ; <br /> RowBox[{RowBox[{a, =, 0.}], ;}] RowBox[{RowBox[{RowBox[{b, =, 2.5}] ... }], }] f[t_] := Sin[t^2] ;  Plot[f[t], {t, a, b}, PlotStyleRed] ;

[Graphics:HTMLFiles/5.9_approxIntegration_4.gif]

Clear[steps,accuracy,left,right,trap,mid,table];

steps = {2,10,50,250};
accuracy = 5;

left[n_] := N[compositeLeftRule[f,a,b,n],accuracy];
right[n_] := N[compositeRightRule[f,a,b,n],accuracy];
trap[n_] := N[compositeTrapezoidRule[f,a,b,n],accuracy];
mid[n_] := N[compositeMidpointRule[f,a,b,n],accuracy];

table = Transpose[{steps,
                   Map[left,steps],
                   Map[right,steps],
                   Map[trap,steps],
                   Map[mid,steps]}];

TableForm[table,
          TableHeadings -> {None,{"n","LEFT","RIGHT","TRAP","MID"}},
             TableSpacing -> {0,3}]

n LEFT RIGHT TRAP MID
2 1.24996 1.20848 1.22922 0.0192431
10 0.461387 0.453092 0.45724 0.416894
50 0.432389 0.43073 0.43156 0.429996
250 0.430725 0.430393 0.430559 0.430497

Let's see what Mathematica gets for the numerical value of this integral.

NIntegrate[f[t],{t,a,b}]

0.430518

Estimating ∫_0^1^(-x^2) x

Hughes-Hallett, Gleason, et al, First Edition, Text Example, page 381

Let's approximate

(* ∫_0^1^(-x^2) x *)

Clear[f, z, a, b, cutPt] <br /> a = 0 ; RowBox[{RowBox[{cutPt, =, RowBox[{1., /, 2^(1/2)}]}],  ...  = 1 ; <br /> f[z_] := Exp[-z^2] ;  Plot[f[z], {z, a, b}, PlotStyleRed] ;

[Graphics:HTMLFiles/5.9_approxIntegration_9.gif]

Clear[n,firstTRAP20,secondMID20,firstMID20,secondTRAP20];

n = 20;

firstTRAP20  = N[compositeTrapezoidRule[f,a,cutPt,n],6];
secondTRAP20 = N[compositeTrapezoidRule[f,cutPt,b,n],6];

firstMID20  = N[compositeMidpointRule[f,a,cutPt,n],6];
secondMID20 = N[compositeMidpointRule[f,cutPt,b,n],6];

Print[firstTRAP20 + secondMID20," <= integral <= ",firstMID20 + secondTRAP20]

RowBox[{0.746734, ,  <= integral <= , , 0.746871}]

Let's see what Mathematica gets for the numerical value of this integral.

NIntegrate[f[z],{z,a,b}]

0.746824

Estimating ∫_0^4 (100 + x^3)^(1/2) x

Hughes-Hallett, Gleason, et al, First Edition, Exercise 7.6.1, page 382

Let's approximate

(* ∫_0^4 (100 + x^3)^(1/2) x *)

Clear[f, x, a, b] <br /> a = 0 ; b = 4 ; <br /> f[x_] := (100 + x^3)^(1/2) ; <br /> Plot[f[x], {x, a, b}, PlotStyleRed] ;

[Graphics:HTMLFiles/5.9_approxIntegration_15.gif]

Clear[steps,accuracy,left,right,trap,mid,table];

steps = {1,2,4};
accuracy = 6;

left[n_] := N[compositeLeftRule[f,a,b,n],accuracy];
right[n_] := N[compositeRightRule[f,a,b,n],accuracy];
trap[n_] := N[compositeTrapezoidRule[f,a,b,n],accuracy];
mid[n_] := N[compositeMidpointRule[f,a,b,n],accuracy];

table = {Map[left,steps],
         Map[right,steps],
         Map[trap,steps],
         Map[mid,steps]};

TableForm[table,
          TableHeadings -> {{"LEFT","RIGHT","TRAP","MID"},
                            {"n=1","n=2","n=4"}},
             TableSpacing -> {0,3}]

n=1 n=2 n=4
LEFT 40. 40.7846 41.7116
RIGHT 51.225 46.3971 44.5179
TRAP 45.6125 43.5909 43.1147
MID 41.5692 42.6386 42.8795

Let's see what Mathematica gets for the numerical value of this integral.

NIntegrate[f[x],{x,a,b}]

42.9581

Estimating ∫_0^1Sin[t^2/2] x

Hughes-Hallett, Gleason, et al, First Edition, Exercise 7.6.5, page 382

Let's approximate

(* ∫_0^1Sin[t^2/2] x *)

Clear[f, t, a, b] <br /> a = 0 ; b = 1 ; <br /> f[t_] := Sin[t^2/2] ;  Plot[f[t], {t, a, b}, PlotStyleRed] ;

[Graphics:HTMLFiles/5.9_approxIntegration_20.gif]

Clear[steps,accuracy,left,right,trap,mid,table];

steps = {10,100};
accuracy = 6;

left[n_] := N[compositeLeftRule[f,a,b,n],accuracy];
right[n_] := N[compositeRightRule[f,a,b,n],accuracy];
trap[n_] := N[compositeTrapezoidRule[f,a,b,n],accuracy];
mid[n_] := N[compositeMidpointRule[f,a,b,n],accuracy];

table = {Map[left,steps],
         Map[right,steps],
         Map[trap,steps],
         Map[mid,steps]};

TableForm[table,
          TableHeadings -> {{"LEFT","RIGHT","TRAP","MID"},
                            {"n=10","n=100"}},
             TableSpacing -> {0,3}]

n=10 n=100
LEFT 0.140474 0.161324
RIGHT 0.188417 0.166118
TRAP 0.164446 0.163721
MID 0.163348 0.16371

Let's see what Mathematica gets for the numerical value of this integral.

NIntegrate[f[t],{t,a,b}]

0.163714

Exercise: Comparison of Several
Numerical Integration Techniques (L, R, T, M)

Exercise: Comparing Several Numerical Integration Techniques

Your task is to approximate the integral, I, of y = Sqrt[x]
for x in the interval [0,1] using various composite methods
for numerical integration. To facilitate using calculators for
these comparisons, we will standardize on using just two
subdivisions of the interval of integration.

Clear[f, x, a, b] ; <br /> f[x_] := x^(1/2)  a = 0 ; b = 1 ; <br /> graph = Plot[f[x], ... #62371;PlotLabely =  Sqrt[x], AxesLabel {"x", "y"}] ;

[Graphics:HTMLFiles/5.9_approxIntegration_23.gif]

Composite Left Formula

Use the 2-fold lefthand-endpoint rule to approximate the integral I.

Illustrate this method by drawing an appropriate figure on the following graph.

Is your result an overestimate or an underestimate of the true value of the integral? Why?

Show[graph];

[Graphics:HTMLFiles/5.9_approxIntegration_24.gif]

Solution

Clear[left,n,deltaX,x1];

n = 2;
deltaX = (b - a)/n //N;
x1 = a + deltaX;

left = deltax (f[a] + f[x1])

RowBox[{0.707107,  , deltax}]

Composite Right Formula

Use the 2-fold righthand-endpoint rule to approximate the integral I.

Illustrate this method by drawing an appropriate figure on the following graph.

Is your result an overestimate or an underestimate of the true value of the integral? Why?

Show[graph];

[Graphics:HTMLFiles/5.9_approxIntegration_26.gif]

Solution

Clear[right];

right = deltaX (f[x1] + f[b])

0.853553

Composite Midpoint Formula

Use the 2-fold midpoint rule to approximate I.

Illustrate this method by drawing an appropriate figure on the following graph.

Is your result an overestimate or an underestimate of the true value of the integral? Why?

Show[graph];

[Graphics:HTMLFiles/5.9_approxIntegration_28.gif]

Solution

Clear[mid];

m1 = a + deltaX/2;
m2 = m1 + deltaX;

mid = deltaX (f[m1] + f[m2])

0.683013

Composite Trapezoid Formula

Use the 2-fold trapezoid rule to approximate I.

Illustrate this method by drawing an appropriate figure on the following graph.

Is your result an overestimate or an underestimate of the true value of the integral? Why?

Show[graph];

[Graphics:HTMLFiles/5.9_approxIntegration_30.gif]

Solution

Clear[trap];

trap = deltaX (f[a] + 2 f[x1] + f[b])

1.20711

Relative Errors

Construct a figure which provides a visual comparison of the actual errors involved
in making the above four approximations. Rank these errors in decreasing order.
There might be a bit of a challenge here in trying to find a picture of the error in the
Midpoint Method. How could we do that?

Show[GraphicsArray[{{graph,graph},
                    {graph,graph}}]];

[Graphics:HTMLFiles/5.9_approxIntegration_32.gif]

Adding Simpson's Rule to the Collection (L, R, T, M, S)

Estimating ∫_1^21/xx

Hughes-Hallett, Gleason, et al, First Edition, Exercise 7.7.2, page 389

Let's estimate

(* ∫_1^21/xx *)

Clear[f,x,a,b,steps,accuracy]

f[x_] := 1/x;

a = 1.;
b = 2.;

steps = 10;
accuracy = 10;

simp10 = N[compositeSimpsonsRule[f,a,b,steps],accuracy];
error  = N[simp10 - Log[2],accuracy];

Print["Simpson's Rule with n = 10 gives  ",simp10]
Print["The true value of the integral is ",Log[2]//N]
Print["              ... so the error is ",simp10 - Log[2]//N]

RowBox[{Simpson's Rule with n = 10 gives  , , 0.693147}]

RowBox[{The true value of the integral is , , 0.693147}]

RowBox[{              ... so the error is , , 1.94105*10^-7}]

Estimating ∫_0^17x^6x

Hughes-Hallett, Gleason, et al, First Edition, Exercise 7.7.11, page 390

Let's estimate

(* ∫_0^17x^6x *)

Clear[f, t, a, b] <br /> a = 0 ; b = 1 ; <br /> f[x_] := 7x^6 ;  Plot[f[x], {x, a, b}, PlotRangeAll, PlotStyleRed] ;

[Graphics:HTMLFiles/5.9_approxIntegration_41.gif]

Clear[steps,accuracy,left,right,trap,mid,simp,table];

steps = {5,10};
accuracy = 6;

left[n_] := N[compositeLeftRule[f,a,b,n],accuracy];
right[n_] := N[compositeRightRule[f,a,b,n],accuracy];
trap[n_] := N[compositeTrapezoidRule[f,a,b,n],accuracy];
mid[n_] := N[compositeMidpointRule[f,a,b,n],accuracy];
simp[n_] := N[compositeMidpointRule[f,a,b,2 n],accuracy];

table = Transpose[{steps,
                   Map[left,steps],
                   Map[right,steps],
                   Map[trap,steps],
                   Map[mid,steps],
                   Map[simp,steps]}];

TableForm[table,
          TableHeadings -> {None,{"n","LEFT","RIGHT","TRAP","MID","SIMP"}},
             TableSpacing -> {0,3}]

n LEFT RIGHT TRAP MID SIMP
5 0.438144 1.83814 1.13814 0.931623 0.982602
10 0.684884 1.38488 1.03488 0.982602 0.995631

Let's see what Mathematica gets for the numerical value of this integral.

NIntegrate[f[t],{t,a,b}]

1.

Estimating ∫_ (-10)^10^(-x^2/2) x

Hughes-Hallett, Gleason, et al, First Edition, Exercise 7.7.11, page 390

Let's estimate

(* ∫_ (-10)^10^(-x^2/2) x *)

Clear[f, t, a, b] <br /> a = -10 ; b = 10 ;  f[x_] := ^(-x^2/2) ; <br /> Plot[f[x], {x, a, b}, PlotRangeAll, PlotStyleRed] ;

[Graphics:HTMLFiles/5.9_approxIntegration_46.gif]

Clear[steps,accuracy,left,right,trap,mid,simp,table];

steps = {5,10,50,100};
accuracy = 6;

left[n_] := N[compositeLeftRule[f,a,b,n],accuracy];
right[n_] := N[compositeRightRule[f,a,b,n],accuracy];
trap[n_] := N[compositeTrapezoidRule[f,a,b,n],accuracy];
mid[n_] := N[compositeMidpointRule[f,a,b,n],accuracy];
simp[n_] := N[compositeMidpointRule[f,a,b,2 n],accuracy];

table = Transpose[{steps,
                   Map[left,steps],
                   Map[right,steps],
                   Map[trap,steps],
                   Map[mid,steps],
                   Map[simp,steps]}];

TableForm[table,
          TableHeadings -> {None,{"n","LEFT","RIGHT","TRAP","MID","SIMP"}},
             TableSpacing -> {0,3}]

n LEFT RIGHT TRAP MID SIMP
5 1.08268 1.08268 1.08268 4.00268 2.47057
10 2.54268 2.54268 2.54268 2.47057 2.50663
50 2.50663 2.50663 2.50663 2.50663 2.50663
100 2.50663 2.50663 2.50663 2.50663 2.50663

Let's see what Mathematica gets for the numerical value of this integral.

NIntegrate[f[t],{t,a,b}]

2.50663


Created by Mathematica  (April 21, 2004)