MODULE 1. 25
Exercises—due 8:00 am Friday 16th August Required
The grade of 1–5 awarded for these exercises is based on the required exercises. If all of the required exercises are completed correctly, a grade of 5 will be obtained.
R1.1 The intensity of a dougnut-shaped (in cross-section) laser beam was mea- sured. The data had the background noise subtracted, was smoothed, and left and right side data (which should be identical) was averaged, giving the intensity as a function of radius:
160
140
120
100
80
60
40
20
0
0 5 10 15 20 25 30
The xy data is:
x=[0 0.1341 0.2693 0.4034 0.5386 0.6727 … 0.8079 0.9421 1.0767 1.2114 1.3460 1.4801 … 1.6153 1.7495 1.8847 2.0199 2.1540 2.2886 … 2.4233 2.5579 2.6921 2.8273 2.9614 3.0966 … 3.2307 3.3659 3.5000 ];
y=[0 0.0310 0.1588 0.3767 0.6452 0.8780 … 0.9719 1.0000 0.9918 0.9329 0.8198 0.7707 … 0.8024 0.7674 0.6876 0.5937 0.5778 0.4755 … 0.3990 0.3733 0.2870 0.2156 0.2239 0.1314 … 0.1180 0.0707 0.0259 ];
Note that the data has been re-scaled for the graph shown; your graph will have different values but should have the same shape.
Provide a simple approximation to this data. The one obtained in the actual real-life case is shown below; this was a piece-wise polynomial approxi- mation. You can use either equation (1.1) or trial-and-error and your own judgment to find the best order of polynomial to use.
26 COSC2500/COSC7500—Semester 2, 2019
160
140
120
100
80
60
40
20
0
−20
0 5 10 15 20 25 30
R1.2 For the polynomial f (x) = x4 − 2×3, calculate the derivative and the 2nd derivative analytically (the first derivative so that you can compare with your numerical derivatives to determine the error, and the second deriva- tive might be useful for understanding your results—see equations (1.9) and (1.11)).
a) Choose different step sizes h and plot a log-log graph of error vs step size to compare the numerical derivatives found using the forward, backward, and central difference formulas with the actual derivative. How accurate are the numerical derivatives at x = 0, 1, 1.5, and 2? At what step sizes do you obtain the minimum errors, and what are the minimum errors (see Sauer figure 5.1)? Do you expect that this would be a problem? What if you were using single-precision floating point? (You can easily try this in Matlab by using the single() function.)
b) Add some additional error to your polynomial. For example, try y = f(x) + 0.001 * randn(size(x));. How does the error in the numer- ical derivative vary with step size h? Try changing the size of the error.
R1.3 Compare the trapezoid rule and Simpson’s rule for integrating cos(x) from x = 0 to x = π/2 and x = π, and sin(x2) + 1 from x = 0 to x = 10. You can use the Matlab function trapz and the Simpson’s method function SIMP42 from Blackboard, or write your own code, or use other code.
a) Plot a log-log graph to compare the accuracy of the trapezoid rule and
Simpson’s rule for different step sizes. 10 2
Note that 0 sin(x ) + 1dx ≈ 10.58367089992962334
b) Compare the accuracy of an adaptive step size integrator (e.g., integral() or quad() in Matlab).
c) How does additional error affect integration? Try it and see, using the same random error term you used before.
d) How much does round-off error affect integration? Since a small value of h means that we need many points, we also increase our require- ments in time and memory. Do the time and/or memory required be- come excessive before we run into any problems with round-off error?
MODULE 1. 27
WARNING: Using too much memory can cause Matlab and/or your computer to crash or become unrespnsive. Don’t increase the number of steps too quickly when trying this (doubling the number of steps should be OK), and make sure you have saved your work.
R1.4 Integrate exp(−x) from x = 0 to x = ∞, by numerically integrating from x = 0 to increasing values of x. Plot a graph of the error vs the value of x you integrate to. Compare results for a fixed step size and a fixed number of steps.
Additional
Attempts at these exercises can earn additional marks, but will not count towards the grade of 1–5 for the exercises. Completing all of these exercises does not mean that 4 marks will be obtained—the marks depend on the quality of the answers. It is possible to earn all 4 marks without completing all of these additional exercises.
A1.1 Howarethespecialfunctionsavailableonascientificcalculatorcalculated?
A1.2 While Matlab is an interpreted language, and therefore loops can be slow, operations on matrices and vectors can be much faster. Matlab code is often written to make use of vector operations instead of loops (“vectorisation”); this can lead to code that is fast, but requires much more memory. Compare the memory and time requirements for numerical integration in Matlab us- ing vector operations, Matlab using loops, and a compiled language such as C using loops.
A1.3 Writeaprogramtocalculatetheincreaseintheviscousdragcoefficientdue to a wall using the approximation given by Chaoui and Feuillebois (2003).
A1.4 ConsiderthealgorithmforthecalculationofWigner3jfunctionsdescribed by Brock (2001). Note that there are a number of special cases. How impor- tant is it to treat these cases separately? A Matlab implementation of Brock’s algorithm, wigner3j.m is available on the course website. (If you compare results with the examples given in Brock (2001), note that the “baseball” calculation is wrong, due to integer overflow.)
A1.5 Estimate the error in a numerical derivative and/or integral from the con- vergence of the solution as the step size is changed. (It’s best to use step sizes that are large enough to avoid round-off error.) Then compare your estimate of the error with actual error determined from comparison with the analytical solution.
A1.6 How would you calculate the integral ∞ 3/2
F(β)= 0 exp(−y −αy)ysin(βy)dy
where α varies from approximately 0.1 to 1, from β = 0.01 to β = 100?
What difficulties do you expect?
Calculate this integral for various combinations of α and β.
28 COSC2500/COSC7500—Semester 2, 2019
A1.7 Using Fourier transforms, the n-th derivative can be approximated by f(n)(x) = F−1 (2πik)n F[f(x)],
where F and F−1 denote the Fourier transform and inverse Fourier trans- form. By replacing F and F−1 with the discrete Fourier transform and its inverse we can write
n 1 N−1 ̄ n N−1 −2πi jk/N 2πi jk/N f(x)=N∑ 2πik ∑f(x)e e ,
k=0 j=0
wherek ̄ =[0,1,2,…N/2−1,N/2,−N/2+1,−N/2+2,…−1].
Implement numerical differentiation using discrete Fourier transforms. The transform and its inverse are available in Matlab: fft and ifft (use OpenCV and Armadillo if you prefer C/C++). When might a Fourier transform approach be useful? For what boundary conditions is the discrete Fourier transform method most accurate?
(Thanks to Isaac Lenton for this question!)
Programming hints
R1.1 It’s useful to plot your fitted curve over a much finer grid of x values rather than just using the original x values. For example, after using polyfit to find your polynomial:
x2 = linspace (0 ,3.5 ,1000); y2 = polyval(P,x2); plot(x,y,’*’,x2,y2);
Now you can see what your fitted curve does between the original data points. (Why is this important?)
R1.2 Sinceyouwanttolookattheeffectofstepsizeonerror,youwanttofindthe numerical derivatives for a many step sizes, over a large range. It’s useful to use logarithmically spaced points to make a vector of step sizes. We have many different ways to make a vector of, e.g., step sizes. The three main ones are:
– Linear spacing, with a specified spacing: h = start_val:spacing:end_val;
– Linear spacing, with a specified number of points:
h = linspace(start_val,end_val,num_points);
– Logarithmic spacing, with a specified number of points:
h = logspace(start_power ,end_power ,num_points);
Note that we don’t specify the start and end values, but the powers of 10.
MODULE 1. 29
In this case, the third one is our simplest choice.
Rather than hard-coding the function you are differentiating in your main script or function, it’s often useful to make a separate m-file for that func- tion, or use an inline function. Very useful when you want to change it, since there’s only one place where you have to make changes! If we wanted a separate m-file to calculate y = x2 exp(x), we could put the following:
function y = f12(x) y = x.^2 .* exp(x); return
in an m-file f12.m (the function name and the m-file name need to be the same). In the first line, we tell Matlab that this is a function (rather than a script), say that there is one output variable (y), and one input parameter (x).
If our function is simple enough to have on one line, we might use an in-line function instead. E.g.,
f12 = @(x) x.^2 .* exp(x);
Here, f12 is a function handle (like a pointer or address). The @ op- erator creates the handle. For more on in-line functions, see Ap- pendix B of Sauer or http://au.mathworks.com/help/matlab/matlab prog/anonymous-functions.html
Next, you need to cycle through the different step sizes. It’s much easier to do this automatically in a loop rather than choosing different step sizes by hand each time. Your code might look something like:
x = 1;
num_points = 1000;
h =
for
%
logspace(-20,-1,num_points);
n = 1:num_points
Forward difference
fd = ( f12(x+h(n)) – f13(x) ) / h(n);
% Exact result
ed = f12exact(x);
% Forward difference error
fde(n) = abs( fd – ed ); end
figure;
loglog(h,fde); title([’x␣=␣’ num2str(x)]); ylabel(’Error’); xlabel(’Step␣size␣h’);
30 COSC2500/COSC7500—Semester 2, 2019 R1.3 To use integral, you need to consider two points:
– Unlike trapz, where the input is vectors of x and y values, the input to integral is a function handle, and start and end x values for the integral. If you have your function in an m-file f13.m, then you need to use the @ operator to produce a function handle for quad:
integral13 = integral(@f13,start_val,end_val);
If you used an inline function for f15, it’s already a function handle, and
integral13 = integral(f13,start_val,end_val);
works.
– You choose how many points trapz and SIMP42 use, while integral
chooses its own number of points. So how can you compare the errors? A simple solution is to find the number of points used by integral. integral will return a 2nd output, the number of function evaluations it used.
CORRECTION: This (2nd output with number of function evalua- tions) doesn’t work with the current version of integral. It does work with quad.