Ma 449: Numerical Applied Mathematics Model Solutions to Homework Assignment 8 Prof. Wickerhauser Exercise 4 of Section 7.1, p.362. The integral of f0*(x-x1)/(x0-x1)+f1*(x-x0)/(x1-x0) is f0*(x-x1)^2/ 2(x0-x1)+f1*(x-x0)^2/2(x1-x0) on [x0,x1] = [f1*(x-x0)^2 - f0*(x-x1)^2] / 2(x1-x0) on [x0,x1] = [f1*(x1-x0)^2 + f0*(x0-x1)^2] / 2(x1-x0) = [f1*(x1-x0)^2 + f0*(x1-x0)^2] / 2(x1-x0) = ( [f1+f0]/2 ) * (x1-x0), which is the trapezoid rule. Exercise 5 of Section 7.1, p.363. Write I for the exact value of the integral and T for the trapezoid rule approximation. Then, for various polynomials p=p(x): p0(x)=1: I(p0) = 1*(x1-x0) T(p0) = ([1+1]/2) *(x1-x0) ==> T(p0)=I(p0) p1(x)=x: I(p1) = [x1^2-x0^2]/2 = ([x1+x0]/2)*(x1-x0) T(p1) = ([x0+x1]/2) *(x1-x0) ==> T(p1)=I(p1) p2(x)=x^2: I(p2) = [x1^3-x0^3]/3 = ([x1^2+x0*x1+x0^2]/3)*(x1-x0) T(p2) = ([x0^2+x1^2]/2) *(x1-x0) ==> T(p2)!=I(p2) Hence the degree of precision of the trapezoid rule is 1. Exercise 6 of Section 7.1, p.363. Write I for the exact value of the integral and S for the Simpson rule approximation. Then, for various polynomials p=p(x): p0(x)=1: I(p0) = 1*(x1-x0) S(p0) = ([1+4*1+1]/6) *(x1-x0) ==> S(p0)=I(p0) p1(x)=x: I(p1) = [x1^2-x0^2]/2 = ([x1+x0]/2)*(x1-x0) S(p1) = ([x0+4*(x0+x1)/2+x1]/6) *(x1-x0) = ([x0+x1]/2) *(x1-x0) ==> S(p1)=I(p1) p2(x)=x^2: I(p2) = [x1^3-x0^3]/3 = ([x1^2+x0*x1+x0^2]/3)*(x1-x0) S(p2) = ([x0^2+4*(x0+x1)^2/4+x1^2]/6) *(x1-x0) = ([2x0^2+2x0*x1+2x1^2]/6) *(x1-x0) = ([x0^2+x0*x1+x1^2]/3) *(x1-x0) ==> S(p2)=I(p2) p3(x)=x^3: I(p3) = [x1^4-x0^4]/4 = ([x1^3+x0*x1^2+x0^2*x1+x0^3]/4)*(x1-x0) S(p3) = ([x0^3+4*(x0+x1)^3/8+x1^3]/6) *(x1-x0) = ([x0^3+x1*x0^2+x1^2*x0+x1^3]/4)*(x1-x0) ==> S(p3)=I(p3) p4(x)=x^4: I(p4) = [x1^5-x0^5]/5 = ([x1^4+x0*x1^3+x0^2*x1^2+x0^3*x1+x0^4]/5)*(x1-x0) S(p4) = ([x0^4+4*(x0+x1)^4/16+x1^4]/6) *(x1-x0) = ([(5/4)x0^4 + lower-order x0 terms]/6)*(x1-x0) ==> S(p4)!=I(p4) Hence the degree of precision of Simpson's rule is 3. Exercise 2(b) of Section 7.2, p.374. C PROGRAM: /* Exercise 2(b) i and ii, chapter 7.2, p.374 */ #include #define PI 3.14159265358979323 typedef double real; typedef real (*integrand)(real); /* Compute the length element for the part (b) function sin(x) */ real length_element(real x) { real c; c = cos(x); /* derivative of sin(x) */ return sqrt(1+c*c); } /* Formula from Th 7.2, Eq 1a, p.364: */ real composite_trapezoid_rule( integrand f, real a, real b, int m ) { real h, sum, x; int k; h = (b-a)/m; sum = (f(a)+f(b))/2; for(k=1; k sqrt( M2 / 12*e ) ~ 2 e 3, so m=2200 should work. Algorithm 2(b) of Section 7.2, p.377. The composite Simpson rule error is bounded by e=M4*(b-a)*h^4/180, where M4 is the maximum value of |f''''(c)| on [a,b], and h=(b-a)/2m. But for f(x)=sqrt(1+cos(x)^2) on [0,pi/4], we have f''''(x) = ( 7290 + 11986 Cos[2 x] + 1128 Cos[4 x] + 45 Cos[6 x] + 9/2 9/2 30 Cos[8 x] + Cos[10 x]) / (2 (3 + Cos[2 x]) ), which is no bigger than 2.2 on [0,pi/4]. [This can be estimated with Matlab's symbolic algebra and plotting routines.] Also, (b-a)=pi/4 and h=pi/8m for this interval. Thus, to satisfy e<5e-12 for 11 decimal places, we must have h^4 < 920*e / M4*pi ==> h < .005, m > pi/8h ~ 77 so m=80 should work. C PROGRAM: /* Algorithm 2 for Exercise 2(b) ii, chapter 7.2, p.377 */ #include typedef double real; typedef real (*integrand)(real); #define PI 3.14159265358979323 /* Compute the length element for the part (b) function sin(x) */ real length_element(real x) { real c; c = cos(x); /* derivative of sin(x) */ return sqrt(1+c*c); } /* Formula from Th 7.3, Eq 4a, p.376: */ real composite_Simpson_rule( integrand f, real a, real b, int m ) { real h, sum, x; int k; h = (b-a)/(2*m); sum = (f(a)+f(b)); for(k=1; k