Download Numerical Computing Problem Set 6 Solutions - Prof. Donald W. Schwendeman and more Assignments Mathematics in PDF only on Docsity!
D. Schwendeman
Numerical Computing Problem Set 6
Due: Thursday, 11/3/
- Consider the integral I =
∫ (^2) h
0
f (x) dx, where h is a constant.
(a) Find the interpolating polynomial f˜ (x) of degree two (or less) fitted to the data (0, f (0)), (h, f (h)) and (2h, f (2h)). Write the polynomial in terms of the Lagrange basis.
(b) Integrate the interpolating polynomial in part (a) to obtain the weights (α 1 , α 2 , α 3 ) in the numerical quadrature formula
I = α 1 f (0) + α 2 f (h) + α 3 f (2h) + R(f )
- Consider the numerical quadrature formula ∫ (^3) h 0
f (x) dx = α 1 f (0) + α 2 f (2h) + R(f )
where h is an arbitrary constant.
(a) Find the weights (α 1 , α 2 ) so that R = 0 when f (x) = 1 and f (x) = x, i.e. when f (x) is a polynomial of degree 1 or less.
(b) The remainder term is R = C hp+2f (p+1)(μ), where μ ∈ [0, 3 h]. Consider monomials f (x) = x^2 , x^3 , etc., to find the positive integer p and the constant C in the remainder term.
- Let I =
∫ (^3)
2
cos(x^2 ) dx. Write a short matlab script to compute approximations for I
using
(a) the composite Simpson rule with 2 subintervals.
(b) the three-point Gaussian quadrature formula.
- The composite midpoint rule for m equally spaced subintervals is
∫ (^) b a
f (x) dx = h
∑^ m i=
f (¯xi) + (b^ −^ a)h
2 24
f ′′(μ), μ ∈ [a, b]
where ¯xi is the midpoint of the ith^ subinterval and h = (b − a)/m. Let f (x) = x/(x + 1), a = − 1 /2 and b = 2. Consider the error term in the composite midpoint rule to determine the number of subintervals needed so that the absolute error in the numerical quadrature is less than 10−^4.
- Let I =
∫ (^) b
a
f (x) dx.
(a) Write a matlab function, mySimpson say, that outputs an approximation for I using the composite Simpson rule with m subintervals. Your function should take input a, b, m and the function f (specified in an M-file).
(b) Let I =
∫ (^2) − 1 exp(−x) sin(5x)^ dx. Find approximations to^ I^ using your matlab function in part (a) with m = 10, 20 and 40. Compute the exact value for I (using Maple if you have to) and use it to find the absolute error in each approximation.
Problem 1
Part (a) We fondly recall Lagrange interpolation from the last homework. The basis functions
are given by
`j (t) =
∏n ∏k=1,k^6 =j^ (t^ −^ tk) n k=1,k 6 =j (tj^ −^ tk)^
Hence 1 (t) = (t(−−hh)()(t−− 22 hh)) , 2 (t) = t((ht)(−−^2 hh)) , and ` 3 (t) = (^) (2t(ht−)(hh)). Then our interpolant is
p(x) = f (0)1 (t) + f (h) 2 (t) + f (2h)` 3 (t) = f^ (0) 2 h^2
(t − h)(t − 2 h) − f^ (h) h^2
t(t − 2 h) + f^ (2h) 2 h^2
t(t − h).
Part (b)
I =
∫ (^2) h
0
f (x)dx =
∫ (^2) h
0
p(x)dx + R(f ) =
f (0) 2 h^2
2 h^3 3 +^
f (h) h^2
4 h^3 3 +^
f (2h) 2 h^2
2 h^3 3 +^ R(f^ )
=⇒ α 1 =
h 3 ,^ α^2 =
4 h 3 ,^ α^3 =^
h 3
Problem 2
Part (a) f (x) = 1 implies
∫ (^3) h 0 f^ (x)dx^ = 3h^ =^ α^1 f^ (0) +^ α^2 f^ (2h) =^ α^1 +^ α^2 ;^ f^ (x) =^ x^ implies ∫ (^3) h 0 f^ (x)dx^ = 9h
(^2) /2 = α 1 f (0) + α 2 f (2h) = 2hα 2. Solving this system for the αi yields α 1 =^3 h 4 ,
α 2 =^94 h.
Part (b) With the αi known, let’s plug in f (x) = x^2 and f (x) = x^3 and find some equations for
R. f (x) = x^2 implies
∫ (^3) h 0 f^ (x)dx^ = 9h
(^3) = α 1 f (0) + α 2 f (2h) + R = 9 h 4 4 h
(^2) + R = 9h (^3) + R, meaning
R = 0. Let’s see what happens with f (x) = x^3. We get
∫ (^3) h 0 f^ (x)dx^ =^
81 h^4 4 =^ α^1 f^ (0)+α^2 f^ (2h)+R^ = 9 h 4 8 h
(^3) + R = 18h (^4) + R. Or, R = ( 81 4 −^ 18)h
(^4) = Chp+2f (p+1)(μ). Matching powers of h, setting
p = 2 yields R = ( 814 − 18)h^4 = Ch^4 f (3)(μ) = Ch^4 · 6, or C =
81 4 −^18 6
=^3
Problem 3
Part (a) Here’s what your code might look like for the general Simpson’s rule.
into the inequality E < 10 −^4 , and solve for m.
We see that f ′′(x) = − (^) (x+1)^23 , whose maximum absolute value over the interval [− 12 , 2] is 16,
occuring at x = − 12. Hence E =
∣ (b−a)
3 24 m^2 f^
′′(μ)
∣ ≤^2.^5
3 24 m^2 16 =^
125 12 m^2 <^10
− (^4) means
m >
3 ≈^322.^75.
Problem 5
Part (a) This is just the code from problem 3.
Part (b) You might run the code using the following function file.
1 function f=myIntegrand2(x) 2 3 % integrand f(x) for mySimpson.m 4 5 f=exp(-x).sin(5x);
And you might use this as a driver script.
1 % Approximate the value of 2 % 3 % I=int(exp(-x)sin(5x), x=-1..2) 4 % 5 % using the composite Simpson rule with m=10, 20 and 40 6 % equally spaced subintervals. 7 8 m=10; 9 S1=mySimpson('myIntegrand2',-1,2,m); 10 11 m=20; 12 S2=mySimpson('myIntegrand2',-1,2,m); 13 14 m=40;
15 S3=mySimpson('myIntegrand2',-1,2,m); 16 17 % Here is the exact value 18 I=0.2732077083; 19 20 % Print results 21 fprintf(1,'I=%0.6e\n',I) 22 fprintf(1,'S(m=10)=%0.6e Error=%0.2e\n',S1,abs(S1-I)) 23 fprintf(1,'S(m=20)=%0.6e Error=%0.2e\n',S2,abs(S2-I)) 24 fprintf(1,'S(m=40)=%0.6e Error=%0.2e\n',S3,abs(S3-I))
The results should be something like this.
I=2.732077e- S(m=10)=2.728739e-001 Error=3.34e- S(m=20)=2.731892e-001 Error=1.85e- S(m=40)=2.732066e-001 Error=1.12e-
