# Numerical Quadrature# # # We test some of our numerical quadrature rules. We first define some# data points.> restart; Digits := 10;> x0 := 0; x1 := evalf(Pi/8); x2 := evalf(Pi/4); x3 := evalf(3*Pi/8); x4> := evalf(Pi/2);> fx0 := evalf(cos(x0)); fx1 := evalf(cos(x1)); fx2 := evalf(cos(x2));> fx3 := evalf(cos(x3));fx4 := evalf(cos(x4));# The trapezoidal rule arises from degree one Lagrange polynomial# interpolation.> ht := x4 - x0;> trap:= (ht/2)*(fx0+fx4);> int(interp([x0,x4],[fx0,fx4],x),x=x0..x4);> traperr := int(cos(x),x=x0..x4) - trap;> # The theoretical bound on the error is given by> ht^3/12;> plot({cos(x),interp([x0,x4],[fx0,fx4],x)},x=x0..x4);# Move on to Simpson's rule> hs := x2-x0;# Simpson's rule arises from degree two Lagrange polynomial# interpolation> simp := (hs/3)*(fx0+4*fx2+fx4);> int(interp([x0,x2,x4],[fx0,fx2,fx4],x),x=x0..x4);> simperr := int(cos(x),x=x0..x4) - simp;# The theoretical bound on the error in Simpson's rule is> hs^5/90;> plot({cos(x),interp([x0,x2,x4],[fx0,fx2,fx4],x)},x=x0..x4);# We can also compute the 5-point closed Newton-Cotes formula using the# points we have available.> h5 := x1 - x0;> ClNC5 := (2*h5/45)*(7*fx0+32*fx1+12*fx2+32*fx3+7*fx4);> int(interp([x0,x1,x2,x3,x4],[fx0,fx1,fx2,fx3,fx4],x),x=x0..x4);> ClNC5err := int(cos(x),x=x0..x4) - ClNC5;# The theoretical bound on the error in this rule is> 8*h5^7/945;> plot({cos(x),interp([x0,x1,x2,x3,x4],[fx0,fx1,fx2,fx3,fx4],x)},x=x0..x> 4);# We turn now to some open Newton-Cotes approximations. First we compute# the midpoint rule, which amounts to the constant function at the# midpoint> hm := x2 - x0;> midp := 2*hm*fx2;> int(interp([x2],[fx2],x),x=x0..x4);> midperr := int(cos(x),x=x0..x4) - midp;# The theoretical bound for the midpoint rule is> hm^3/3;> plot({cos(x),interp([x2],[fx2],x)},x=x0..x4);# We can also compute the n=2 open Newton-Cotes formula (note we have to# shift the subscripts by one to match our labeling.> h2 := x1-x0;> OpNC2 := (4*h2/3)*(2*fx1-fx2+2*fx3);> int(interp([x1,x2,x3],[fx1,fx2,fx3],x),x=x0..x4);> OpNC2err := int(cos(x),x=x0..x4) - OpNC2;# The theoeretical bound for this rule is> 14*h2^5/45;> plot({cos(x),interp([x1,x2,x3],[fx1,fx2,fx3],x)},x=x0..x4);# Now we turn to the composite rules on our interval> with(student);# First, the composite trapezoid rule:> ct := evalf(trapezoid(cos(x),x=x0..x4,4));> hct := x1-x0;> hcterr := int(cos(x),x=x0..x4) - ct;> (x4-x0)*hct^2/12;# Next, the composite midpoint rule:> hmt := x1-x0;> mt := 2*hmt*(fx1+fx3);> hmterr := int(cos(x),x=x0..x4) - mt;> (x4-x0)*hmt^2/6;# Finally, the composite Simpson's rule> st := evalf(simpson(cos(x),x=x0..x4,4));> hst := x1-x0;> hsterr := int(cos(x),x=x0..x4) - st;> (x4-x0)*hst^4/180;# Starter values for Romberg integration> evalf(trapezoid(cos(x),x=x0..x4,1));> evalf(trapezoid(cos(x),x=x0..x4,2));> evalf(trapezoid(cos(x),x=x0..x4,4));> evalf(trapezoid(cos(x),x=x0..x4,8));> evalf(trapezoid(cos(x),x=x0..x4,16));>