MATH 132 -- Calculus for Physical and Life Science 2

Solutions for Lab Day 1 -- "In Search of a Better Numerical

Integration Method"

A)  1)  First integral:   int(exp(-x^2/10),x = 0 .. 2)

>    with(student);

[D, Diff, Doubleint, Int, Limit, Lineint, Product, Sum, Tripleint, changevar, completesquare, distance, equate, integrand, intercept, intparts, leftbox, leftsum, makeproc, middlebox, middlesum, midpoin...
[D, Diff, Doubleint, Int, Limit, Lineint, Product, Sum, Tripleint, changevar, completesquare, distance, equate, integrand, intercept, intparts, leftbox, leftsum, makeproc, middlebox, middlesum, midpoin...

>    f1:=exp(-x^2/10);

f1 := exp(-1/10*x^2)

>    exact1:=evalf(Int(f1,x=0..2));

exact1 := 1.762508070

>    n:='n';

n := 'n'

>    with(Spread):

>    Table1:=CreateSpreadsheet():

n LEFT MID RIGHT TRAP LEFT*Error MID*Error RIGHT*Error TRAP*Error
5 1.824864814 1.764299271 1.692992833 1.758928823 -.62356744e-1 -.1791201e-2 .69515237e-1 .3579247e-2
10 1.794582043 1.762955179 1.728646052 1.761614048 -.32073973e-1 -.447109e-3 .33862018e-1 .894022e-3
20 1.778768611 1.762619804 1.745800616 1.762284613 -.16260541e-1 -.111734e-3 .16707454e-1 .223457e-3
40 1.770694208 1.762536001 1.754210210 1.762452209 -.8186138e-2 -.27931e-4 .8297860e-2 .55861e-4
80 1.766615104 1.762515052 1.758373106 1.762494105 -.4107034e-2 -.6982e-5 .4134964e-2 .13965e-4
160 1.764565079 1.762509815 1.760444079 1.762504579 -.2057009e-2 -.1745e-5 .2063991e-2 .3491e-5

Table1 := SpreadSheet007

>    ns:=[5,10,20,40,80,160];

ns := [5, 10, 20, 40, 80, 160]

>    for i to 6 do
  SetCellFormula(Table1,i+1,1,ns[i]);
  SetCellFormula(Table1,i+1,2,evalf(leftsum(f,x=0..2,ns[i])));
  SetCellFormula(Table1,i+1,3,evalf(middlesum(f,x=0..2,ns[i])));
  SetCellFormula(Table1,i+1,4,evalf(rightsum(f,x=0..2,ns[i])));
  SetCellFormula(Table1,i+1,5,evalf(trapezoid(f,x=0..2,ns[i])));
  SetCellFormula(Table1,i+1,6,exact1-evalf(leftsum(f,x=0..2,ns[i])));
  SetCellFormula(Table1,i+1,7,exact1-evalf(middlesum(f,x=0..2,ns[i])));
  SetCellFormula(Table1,i+1,8,exact1-evalf(rightsum(f,x=0..2,ns[i])));
  SetCellFormula(Table1,i+1,9,exact1-evalf(trapezoid(f,x=0..2,ns[i])));

>    end do;

>    EvaluateSpreadSheet(Table1);

EvaluateSpreadSheet(SpreadSheet007)

>   

  2)  Second integral:   int(sin(x)/x,x = 2 .. 5)

>    f2:=sin(x)/x;

f2 := sin(x)/x

>    Table2:=CreateSpreadsheet();

n LEFT MID RIGHT TRAP LEFT*Error MID*Error RIGHT*Error TRAP*Error
5 .1544328713 -.6350025936e-1 -.2334272696 -.3949719918e-1 -.2099146032 .801852750e-2 .1779455377 -.1598453268e-1
10 .4546630599e-1 -.5747486457e-1 -.1484637645 -.5149872926e-1 -.1009480378 .199313271e-2 .9298203264e-1 -.398300260e-2
20 -.6004279299e-2 -.5597930114e-1 -.1029693146 -.5448679694e-1 -.4947745256e-1 .49756928e-3 .4748758274e-1 -.99493492e-3
40 -.3099179022e-1 -.5560607963e-1 -.7947430785e-1 -.5523304903e-1 -.2448994164e-1 .12434777e-3 .2399257599e-1 -.24868283e-3
80 -.4329893494e-1 -.5551281600e-1 -.6754019374e-1 -.5541956433e-1 -.1218279692e-1 .3108414e-4 .1205846188e-1 -.6216753e-4
160 -.5548950272e-1 -.5548950272e-1 -.6152650489e-1 -.5546619018e-1 -.607585639e-2 .777086e-5 .604477303e-2 -.1554168e-4










-.5548173186e-1






Table2 := SpreadSheet008

>   

>    ns:=[5,10,20,40,80,160];

ns := [5, 10, 20, 40, 80, 160]

>    exact2:=evalf(Int(f2,x=2..5));

exact2 := -.5548173186e-1

>    for i to 6 do
  SetCellFormula(Table2,i+1,1,ns[i]);
  SetCellFormula(Table2,i+1,2,evalf(leftsum(f2,x=2..5,ns[i])));
  SetCellFormula(Table2,i+1,3,evalf(middlesum(f2,x=2..5,ns[i])));
  SetCellFormula(Table2,i+1,4,evalf(rightsum(f2,x=2..5,ns[i])));
  SetCellFormula(Table2,i+1,5,evalf(trapezoid(f2,x=2..5,ns[i])));
  SetCellFormula(Table2,i+1,6,exact2-evalf(leftsum(f2,x=2..5,ns[i])));
  SetCellFormula(Table2,i+1,7,exact2-evalf(middlesum(f2,x=2..5,ns[i])));
  SetCellFormula(Table2,i+1,8,exact2-evalf(rightsum(f2,x=2..5,ns[i])));
  SetCellFormula(Table2,i+1,9,exact2-evalf(trapezoid(f2,x=2..5,ns[i])));

>    end do;

>    EvaluateSpreadSheet(Table2);

EvaluateSpreadSheet(SpreadSheet008)

3)  Third Integral:   int(sqrt(1+x^4),x = 0 .. 1)

Table1 := SpreadSheet004

>    ns:=[5,10,20,40,80,160];

ns := [5, 10, 20, 40, 80, 160]

>    f3:=sqrt(1+x^4);

f3 := (1+x^4)^(1/2)

>    exact3:=evalf(Int(f3,x=0..1));

exact3 := 1.089429413

>    Table3:=CreateSpreadsheet();
n LEFT MID RIGHT TRAP LEFT*Error MID*Error RIGHT*Error TRAP*Error
5 1.052722196 1.087072300 1.135564908 1.094143552 .36707217e-1 .2357113e-2 -.46135495e-1 -.4714139e-2
10 1.069897248 1.088840156 1.111318604 1.090607926 .19532165e-1 .589257e-3 -.21889191e-1 -.1178513e-2
20 1.079368702 1.089282100 1.100079380 1.089724041 .10060711e-1 .147313e-3 -.10649967e-1 -.294628e-3
40 1.084325401 1.089392585 1.094680740 1.089503071 .5104012e-2 .36828e-4 -.5251327e-2 -.73658e-4
80 1.086858993 1.089420206 1.092036662 1.089447828 .2570420e-2 .9207e-5 -.2607249e-2 -.18415e-4
160 1.088139599 1.089427111 1.090728434 1.089434016 .1289814e-2 .2302e-5 -.1299021e-2 -.4603e-5

Table3 := SpreadSheet009

>    for i to 6 do
  SetCellFormula(Table3,i+1,1,ns[i]);
  SetCellFormula(Table3,i+1,2,evalf(leftsum(f3,x=0..1,ns[i])));
  SetCellFormula(Table3,i+1,3,evalf(middlesum(f3,x=0..1,ns[i])));
  SetCellFormula(Table3,i+1,4,evalf(rightsum(f3,x=0..1,ns[i])));
  SetCellFormula(Table3,i+1,5,evalf(trapezoid(f3,x=0..1,ns[i])));
  SetCellFormula(Table3,i+1,6,exact3-evalf(leftsum(f3,x=0..1,ns[i])));
  SetCellFormula(Table3,i+1,7,exact3-evalf(middlesum(f3,x=0..1,ns[i])));
  SetCellFormula(Table3,i+1,8,exact3-evalf(rightsum(f3,x=0..1,ns[i])));
  SetCellFormula(Table3,i+1,9,exact3-evalf(trapezoid(f3,x=0..1,ns[i])));

>    end do;

>    EvaluateSpreadSheet(Table3);

EvaluateSpreadSheet(SpreadSheet009)

B)  1)  You might notice that the error decreases when we double n in every case.  Moreover, the

           size of the error seems to decrease roughly by a factor of 2 when we double n for the

           LEFT and RIGHT methods, and roughly by a factor of 4 when we double   n   for the

           MID and TRAP methods.

     2)  When we fix   n   the LEFT and RIGHT errors are usually close to each other in magnitude, and

          also usually quite a bit larger than the  MID and TRAP errors.  Also, when we compare

          the MID and TRAP errors, we see that they have opposite signs, and the TRAP error is

           roughly twice as large in magnitude.

     3)  The MID method gives an overestimate (negative error) if the graph is concave down

          and an underestimate (positive error) if the graph is concave up.  We have the first

          case in integral 1, and the second for integrals 2 and 3

C)  1)  Using Simpson's Rule

.

>    Table4:=CreateSpreadsheet()
n Integral1 Integral1*Error Integral2 Integral2Error Integral3 Integral3Error
10 1.762509122 -.1052e-5 -.5549923929e-1 .1750743e-4 1.089429384 .29e-7
20 1.762508136 -.66e-7 -.5548281947e-1 .108761e-5 1.089429413 0.
40 1.762508073 -.3e-8 -.5548179972e-1 .6786e-7 1.089429413 0.
80 1.762508071 -.1e-8 -.5548173611e-1 .425e-8 1.089429413 0.
160 1.762508070 0. -.5548173210e-1 .24e-9 1.089429413 0.














Table4 := SpreadSheet010

>    for i from 2 to 6 do
  SetCellFormula(Table4,i,1,ns[i]);
  SetCellFormula(Table4,i,2,evalf(simpson(f1,x=0..2,ns[i])));
  SetCellFormula(Table4,i,4,evalf(simpson(f2,x=2..5,ns[i])));
  SetCellFormula(Table4,i,6,evalf(simpson(f3,x=0..1,ns[i])));
 
  SetCellFormula(Table4,i,3,exact1-evalf(simpson(f1,x=0..2,ns[i])));
  SetCellFormula(Table4,i,5,exact2-evalf(simpson(f2,x=2..5,ns[i])));
  SetCellFormula(Table4,i,7,exact3-evalf(simpson(f3,x=0..1,ns[i])));
  end do;

        Note that the errors here are all significantly smaller than the corresponding errors for

        LEFT, RIGHT, and even MID and TRAP.  

  2)   The reason that the SIMPSON method is significantly more accurate is the relationship

        between the MID and TRAP errors we noted in B 3.  Recall that we said that the TRAP

        error was roughly twice as large in magnitude, and of the opposite sign than the

        MID error.   This means when we form the weighted average:

           SIMPSON(n) = (TRAP(n)+2*MID(n))/3

        the errors are almost cancelling out.  (They would exactly cancel out and yield an error

         of zero if the TRAP error was exactly twice as large as the MID error in magnitude).