353 Math Engineering Mathematics III 05F R. J. Braun http://www.math.udel.edu/braun/M353/M353.html Assignment 6, due 9/28/05, 2:30pm 1. Download the TestMyExp script files and the related function files. Run the script files, record the results in a diary file and comment on the relative speed of the scalar and vector versions. Also comment on the ways used to time the approaches. Note: You will likely have to adjust the value of nRepeat in the script files to get a good value of the relative time required that doesn't take too long to run on your computer. In your comments, tell me about the kind of computer you are using: cpu, cache memory, memory, operating system and so forth. 2. Modify the trapezoidal rule script file on the course web page (given to you by Dr. Monk in my absense as well) so that it uses a function file to evaluate the function at every point. (Create a separate function file for this purpose.) Use the code to do the following integrals using n = 16 intervals with the trapezoidal rule. (a). I1 = (b). I2 = (c). I3 =
1 x e dx. 0 3 ln(x)dx. 1 cos(x)dx. -
Compute the error in each case. Solutions, hints and answers Problem 1. Timing Taylor Polynomials for Approximating Exponentials. The following diary file shows some timing results from using TestMyExpRB1.m and TestMyExpRB2.m. >> >> >> >> >> >> % Testing the scalar implementation in MyExp1.m vs the vector implemenation % in MyExp2.m using two different timing commands. % TestMyExpRB1.m uses tic and toc to find wall clock time % TestMyExpRB2.m uses cputime to find actual cpu time % type MyExp1.m
function y = MyExp1(x) % y = MyExp1(x) % x is a column n-vector and y is a column n-vector with the property that % y(i) is a Taylor approximation to exp(x(i)) for i=1:n. nTerms = 17; n = length(x); y = ones(n,1); for i=1:n y(i) = MyExpF(x(i),nTerms); end >> type MyExp2.m function y = MyExp2(x) % y = MyExp2(x) % x is an n-vector and y is an n-vector with with the same shape 1
% and the property that y(i) is a Taylor approximation to % exp(x(i)) for i=1:n. y = ones(size(x)); nTerms = 17; term = ones(size(x)); for k=1:nTerms term = x.*term/k; y = y + term; end >> TestMyExpRB1 % using nRepeat = 50 T1 = MyExp1 benchmark. T2 = MyExp2 benchmark. Length(x) T2/T1 ----------------------1000 0.147430 1100 0.131516 1200 0.138785 1300 0.126010 1400 0.124822 1500 0.139547 1600 0.125225 1700 0.132961 1800 0.132659 1900 0.133946 2000 0.128568 >> TestMyExpRB2 % using nRepeat = 50 T1 = MyExp1 benchmark. T2 = MyExp2 benchmark. Length(x) T2/T1 ----------------------1000 0.125000 1100 0.166667 1200 0.105263 1300 0.150000 1400 0.136364 1500 0.125000 1600 0.160000 1700 0.153846 1800 0.142857 1900 0.133333 2000 0.129032 >> TestMyExpRB2 % using nRepeat = 100 T1 = MyExp1 benchmark. T2 = MyExp2 benchmark. Length(x) T2/T1 ----------------------1000 0.156250 1100 0.147059 1200 0.135135
2
1300 0.150000 1400 0.136364 1500 0.127660 1600 0.140000 1700 0.129630 1800 0.142857 1900 0.135593 2000 0.145161 >> TestMyExpRB2 % using nRepeat = 10 T1 = MyExp1 benchmark. T2 = MyExp2 benchmark. Length(x) T2/T1 ----------------------1000 0.000000 1100 0.333333 1200 0.000000 1300 0.000000 1400 0.200000 1500 0.000000 1600 0.000000 1700 0.200000 1800 0.166667 1900 0.142857 2000 0.166667 >> % using a large value of nRepeat eliminates bogus timings and >> % helps minimize fluctuations in the ratio of times >> diary off The value of nRepeat must be large enough to get reasonable timings for the two methods; a ratio of zero is clearly bogus, and with enough repetitions, the values of T2/T1 stabilize to sensible values. In all cases, the vectorized method of MyExp2.m is faster than MyExp1.m because the ratio is never larger than one. The two timing methods give similar results. The commands tic and toc give a wall-clock time, that is, what you might get if you were very good with a stop watch. The cputime command will count only the time that the cpu spends on the Matlab task at hand, not other processes that the cpu may be doing during the calculation. Solutions, hints and answers Problem 2. First Approach. The following diary file shows a script Trap1.m file that uses three function files, one for each part of the problem. The user inputs the function name (without the .m) and the limits of integration; the script selects the correct exact answer to use with an if-else-end structure. The results show that the first two integrals are approximated well, around 10-4 error or so, while the third one seems particularly good, being about the size of the unit round with only 16 points. We'll discuss this later; something special happens in that last integral due to symmetry. >> type Trap1.m % Script Trap1.m uses the composite trapezoidal rule % to approximate the integral of function f % % Input: f - function to integrate 3
% a - lower limit of integration % b - upper limit of integration % f = input('Type in the function name within single quotes: ') a = input('Type in the lower limit: ') b = input('Type in the upper limit: ') n = 16; % fixed at n = 16 intervals for this problem h = (b-a)/n; % n intervals means n+1 points; use that in linspace x = linspace(a,b,n+1); F = feval(f,x); % do integration now; uses vectorized calculation S = (F(1)+F(n+1))/2 + sum(F(2:n)); ApproxI = S*h % choose between three exact answers if (f == 'Hw7F1') Exact = exp(b) - exp(1); % first problem elseif (f == 'Hw7F2') Exact = b*(log(b) - 1) - a*(log(a) - 1); % second problem elseif (f == 'Hw7F3') Exact = sin(b) - sin(a); % third problem end ErrInApprox = abs(Exact - ApproxI) >> type Hw7F1.m function y = Hw7F1(x) % Function for Hw 7, numerical integration problem % First function y = exp(x); >> Trap1 Type in the function name within single quotes: 'Hw7F1' f = Hw7F1 Type in the lower limit: 0 a = 0 Type in the upper limit: 1 b = 1 ApproxI = 1.7188 ErrInApprox = 5.5930e-04 >> type Hw7F2.m function y = Hw7F2(x) % Function for Hw 7, numerical integration problem % Second function y = log(x); >> Trap1 Type in the function name within single quotes: 'Hw7F2'
4
f = Hw7F2 Type in the lower limit: 1 a = 1 Type in the upper limit: 3 b = 3 ApproxI = 1.2950 ErrInApprox = 8.6741e-04 >> type Hw7F3.m function y = Hw7F3(x) % Function for Hw 7, numerical integration problem % Third function y = cos(x); >> Trap1 Type in the function name within single quotes: 'Hw7F3' f = Hw7F3 Type in the lower limit: -pi a = -3.1416 Type in the upper limit: pi b = 3.1416 ApproxI = 1.7439e-16 ErrInApprox = 7.0536e-17 >> diary off Second Approach. This time, we have function file that does the trapezoidal rule for us; the function for the integ...