Download Numerical Integration using Trapezoidal and Simpson's 1/3 Rule in Matlab - Prof. Aldo Ferr and more Assignments Computer Science in PDF only on Docsity!
ME2016, Fall 2008, Dr. Ferri Homework Assignment 11 Due Wednesday, November 26.
Problem 1 (4 points)
Consider the integral = ∫ = −
3 0
I sin(x)dx 1 cos( 3 ) = 1.
(a) Estimate the value of the integral using trapezoidal integration and a step size of h = Δx = 0.5 (n = 6).
Compare the true truncation error E (^) t = (I (^) true - I (^) approx ) with the approximate relation given in equation (21.13) of the text using an average curvature of f ′′^ = -1/2.
f n
b a E (^) a ′′
3
12
(21.13)
(b) Repeat the calculations of part(a) for h = 0.25 (n = 12). By what factor does the error decrease compared to your answer for part (a)? How does that compare with the theory as given in Eq. (21.13)?
(c) Now, estimate the value of the integral using Simpson's 1/3 rule and h = 0.25. Give the % error relative to the exact answer. Also, comparing the error for the trapezoidal integration with h = 0. (answer from part b) and the present error, by what factor has the error decreased? If the Simpson- 1/3 rule results in an improvement over the trapezoidal integration with the same h , does the improvement match the theoretical prediction indicated in a comparison of (21.13) to (21.19)? (You may use f (^4 ) = +1/2).
( 4 ) 4
5
180
f n
b a E (^) a
Note: Assume that the x-spacing is given by h=(b-a)/n, and that the vector of function values is contained in the 1x(n+1) vector f. In keeping with the Matlab indices that must start at 1 (not 0), the formulas for trapezoidal and Simpson’s-1/3 rule must be modified as follows:
Trapezoidal:
∫ =^ + ∑ +
=
n
i
i n
b a f x f x f x
h f x dx 2
( ) 2 ( 1 ) 2 ( ) ( 1 ) »I= h*(sum(f)-f(1)/2–f(n+1)/2)
Simpson’s-1/3:
−
= =
∫ (^ ) 3 ( )^4 ∑ ( )^2 ∑ ( ) (^1 )
1
2 , 4 , 6 3 , 5 , 7
1 n
n
i
i
n
i
i
b a f x f x f x f x
h f x dx
Matlab: s=f(1); for j=2:2:(n-2); s = s + 4f(j)+2f(j+1); end; s = s + 4f(n)+f(n+1); I = sh/
Solution
(a) The integral can be estimated using the following commands:
a=0; b=3; n=6; h=(b-a)/n; x=(0:n)h; f=sin(x); Sum = f(1)/2; for j=2:n; Sum = Sum + f(j); end Sum = Sum + f(n+1)/2; I05 = Sumh
Equivalently, the last 6 lines could be replaced with the command
I05=(sum(f) - f(1)/2 - f(n+1)/2)*h
Or, one could simply use the Matlab trapz command:
I05 = trapz(x,f)
The result (from any of the three approaches) is I05 = 1.94836054248786, with a truncation error of E (^) t = 0.041632, and a percent error of et = 2.092 %. For comparison, the approximate truncation error is given by equation (21.13) of the text
f n
(b a) E (^) a ′′
3
12
(21.13)
Using f ′′^ = -1/2, we obtain Ea = 0.031250, which compares favorably with E (^) t.
(b) The result for h=0.25 can be obtained by replacing n=6 with n=12 in the code above. The result is I025 = 1.97961713985584, with a truncation error of Et = 0.010375, and a percent error of et = 0.521%. Judging from equation (21.13), doubling n should cause the truncation error to fall by a factor of 4, which is very close to what was found. (The actual truncation error fell by a factor of 4.013.)
(c) Simpson's 1/3 rule can be implemented by the following Matlab-translation of the pseudocode given in Figure 21.14:
a=0; b=3; n=12; h=(b-a)/n; x=(0:n)*h; f=sin(x);
Sum=f(1); for j=2:2:(n-2); Sum = Sum + 4f(j) + 2f(j+1); end Sum = Sum + 4f(n) + f(n+1); I025_simp13 = Sumh/
Note that the "Do-Loop" in the pseudocode of Figure 21.14 ranges from i=1 to i = n-2 in steps of 2; in this case, however, f(1) is f(x 0 ), f(2) = f(x 1 ), etc.
Another way of coding up the Simpson 1/3 rule is as follows:
(a) Calculate the approximate value of I from Problem 1 using trapezoidal integration, and using n = 1, 2, 4, 8, 16, and 32 intervals. Calculate the “true error,” Et = I_exact – In, for each value of n; for I_exact, use I_exact = 1-cos(3) (b) Use expression (2.5) to combine your answers from part (a) to obtain improved estimates of the integral. Does the trend of the error in the new estimates decrease like 1/n 4 as do those from Simpson’s 1/3 rule?
Solution (a) The function, as well as a 4-segment trapezoidal integration approximation is shown in the figure below.
(^00) 0.5 1 1.5 2 2.5 3
1
x
f
Trapezoidal Integration, n= 4 Exact integral = 1. Trapezoidal int. = 1. Percent error = 4.
The trend in the error as a function of n is calculated using the following Matlab code:
Ett = []; It = []; N = [1 2 4 8 16 32]'; for j = 1:6; n = N(j); h = 3/n; xh = 0:h:3; fh = sin(xh); I = h(sum(fh(1:(n+1))) - 0.5fh(1) - 0.5*fh(n+1)); It = [It I]; Ett = [Ett (I_exact - I)]; end
plot(N,Ett,'o-') xlabel('Number of intervals, n') ylabel('Truncation error, Et') title('Trapezoidal Integration') grid
The following plot is produced:
(^00 5 10 15 20 25 30 )
1
Number of intervals, n
Truncation error, Et
Trapezoidal Integration
In tabular form, the results are given by
» disp(' n I Et'); » disp(' '); » disp([N It(:) Ett(:)]);
n I Et
1.00000000000000 0.21168001208980 1. 2.00000000000000 1.60208248595098 0. 4.00000000000000 1.89582521065893 0. 8.00000000000000 1.96661743160631 0. 16.00000000000000 1.98415902154761 0. 32.00000000000000 1.98853476901713 0.
Both the figure and the table show that the error falls off rapidly with n. To see if the error falls off like 1/n 2 , we can calculate the ratio of the error in the integral for some particular value of h (or n) to the error in the integral with h/2 (or 2*n):
» Et1_over_Et2 = Ett(1)/Ett(2) %(error with n=1)/(error with n=2) » Et2_over_Et4 = Ett(2)/Ett(3) %(error with n=2)/(error with n=4) » Et4_over_Et8 = Ett(3)/Ett(4) » Et8_over_Et16 = Ett(4)/Ett(5) » Et16_over_Et32 = Ett(5)/Ett(6)
Gives the results :
Et1_over_Et2 = 4. Et2_over_Et4 = 4. Et4_over_Et8 = 4. Et8_over_Et16 = 4. Et16_over_Et32 = 4.
It is seen that doubling n reduces Et by a factor of 4 (approximately).
(b) The following lines of code are used to generate the improved estimates of the integral:
for j = 1:
(b) Use trapezoidal integration with h= σ/20 to determine the probability that a student will receive an A, P (^) A. Since the upper limit of integration is infinity, integrate first over one interval of width h (from ,
x = x + 1.5 σ to x = x + 1.5 σ+ h ) using a single application of trapezoidal integration. Then,
add the contributions from each subsequent h-interval one-by-one. Continue adding intervals until the approximate error e (^) a =100*abs(Inew - I (^) old)/I (^) new < 0.1 %. Give the value of the probability of receiving an A, and give the maximum value of x that was required for the integral to converge.
Solution One way to accomplish this is to create a Matlab function (separate m-file) and then turn it into an anonymous function. I wrote the Matlab function gauss_dist.m that takes a vector of input values, x, along with the mean and standard-deviation values, and generates a vector of output values corresponding to a Gaussian probability density function:
% Function gauss_dist.m, function takes a vector input x and % gives a vector of values corresponding to the probability % density function having mean xbar and standard deviation % sigma % function P = gauss_dist(x, xbar, sigma); % xbar = 70; sigma = 9; q = -0.5((x-xbar)/sigma).^2; P = exp(q)/sqrt(2pi)/sigma;
(a) To evaluate the probability that a student will obtain a B, type the following Matlab commands:
xbar = 70; sigma = 9; func = @(x) gauss_dist(x,xbar,sigma); aB = xbar + 0.4sigma bB = xbar + 1.5sigma Prob_of_B = quad(func,aB,bB)
This yields Prob_of_B = 0.2778, or, there is a 27.78% chance of a student obtaining a B in the class. (The B range is found to be the range 73.60 to 83.50)
To evaluate the probability that a student will obtain a C, type the following Matlab commands:
aC = xbar - 0.75sigma bC = xbar + 0.4sigma Prob_of_C = quad(func,aC,bC)
This yields Prob_of_C = 0.4288, or, there is a 42.88% chance of a student obtaining a C in the class. (The C range is found to be the range 63.25 to 73.60)
To evaluate the probability that a student will obtain a D, type the following Matlab commands:
aD = xbar - 2.2sigma bD = xbar - 0.75sigma Prob_of_D = quad(func,aD,bD)
This yields Prob_of_D = 0.2127, or, there is a 21.27% chance of a student obtaining a D in the class. (The D range is found to be the range 50.20 to 63.25)
(b) The Matlab code for implementing the trapezoidal-integration of the indefinite integral of
p(x) from x + 1. 5 σ = 83.50 to infinity is given below:
ea = 100; es = 0.1; Int_old = 0; h = sigma/20; aA = xbar + 1.5sigma; a = aA; b = a + h; fa = func(a); while abs(ea) > es; b = a + h fb = func(b); Int_new = Int_old + 0.5h(fa + fb); a = b; fa = fb; ea = 100(Int_new - Int_old)/Int_new; Int_old = Int_new; end Prob_of_A = Int_new
The result is Prob_of_A = 0.0665666, or 6.66% (The algorithm required 39 steps before ea < es.) The result compares favorably with the quad command using a sufficiently high upper integration limit:
aA = xbar + 1.5sigma bA = xbar + 5sigma Prob_of_A = quad(func,aA,bA)
which yields a result of Prob_of_A = 0.0668067 (The upper limit is bA = 115.0)
Knowing that the probabilities must add up to 1, the probability of failure is found to be Prob_of_F = 0.014146. Thus, there is only a 1.41% chance of failing.
The probability density function together with the letter-grade cutoff limits is shown in the figure below:
