关键词 > COSC2500/COSC7500
COSC2500/COSC7500—Semester 2, 2022
发布时间：2022-08-17
Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: daixieit
COSC2500/COSC7500—Semester 2, 2022
Exercises
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:
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 ind the best order of polynomial to use.
MODULE 1. COMPUTATIONAL ERROR 27
160
140
120
100
80
60
40
20
0
−20
0 5 10 15 20 25 30
Checklist
Did you it a suitable function or functions to the data?
Did you include a graph comparing your it to the data?
Did you explain how you chose the best it?
Did you include a graph or it measure values illustrating your choice? Did you include code where appropriate?
Did you reference the source of any code you didn’t write?
R1.2 For the polynomial f (x) = x4 − 2x , calculate the derivative and the 2nd3
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 igure 5.1)? Do you expect that this would be a problem? What if you were using single-precision loating 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.
Checklist
Did you differentiate the function (including with added error)? Did you explain how you did it?
Did you plot graphs showing the error vs step size at each point? Are your graphs and axes properly labelled?
Did you discuss whether the results were what you expected? Did you discuss any unusual results?
Did you include code where appropriate?
Did you reference the source of any code you didn’t write?
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.
Hint: Note that only certain step sizes are allowed. Rather than choosing the step size, choose the number of points, and calculate the step sizefrom that.
Note that \010 sin(x2 ) + 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? 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.
Checklist
Did you integrate the functions using appropriate integrators? Did you explain how you did it?
Did you plot graphs showing the error vs step size for each function? Are your graphs and axes properly labelled?
Did you discuss whether the results were what you expected? Did you discuss any unusual results?
Did you discuss the effect of additional error?
Could you see any effect due to round-off error? Did you discuss this? Did you include code where appropriate?
Did you reference the source of any code you didn’t write?
R1.4 Integrate exp(−x) from x = 0 to x = ∞, by numerically integrating from
Did you integrate the function to increasing values of x?
Did you explain how you did it?
Did you plot graphs showing the error, for ixed step size and step number? Are your graphs and axes properly labelled?
Did you discuss whether the results were what you expected?
Did you discuss any unusual results?
Did you include code where appropriate?
Did you reference the source of any code you didn’t write?
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 How are the special functions available on a scientiic calculator calculated?
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 Write a program to calculate the increase in the viscous drag coeficient due to a wall using the approximation given by Chaoui and Feuillebois (2003).
A1.4 Consider the algorithm for the calculation of Wigner 3j functions described
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 overlow.)
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
F(β) = \0 ∞ exp(−y3/2 − αy)y sin(βy)dy
where α varies from approximately 0.1 to 1, from β = 0.01 to β = 100? What dificulties do you expect?
Calculate this integral for various combinations of α and β.
A1.7 Using Fourier transforms, the n-th derivative can be approximated by f(n)(x) = F − 1 [(2πik)n F [f(x)]J ,
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
fn (x) = 
( (2πi
) n 
f(x)e−2πi jk/N) e2πi jk/N, where 
= [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 itted curve over a much iner grid of x values rather
than just using the original x values. For example, after using polyfit to ind your polynomial:
Now you can see what your itted curve does between the original data points. (Why is this important?)
R1.2 Since you want to look at the effect of step size on error, you want to ind the
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 speciied spacing:
– Linear spacing, with a speciied number of points:
– Logarithmic spacing, with a speciied 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.
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-ile 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-ile to calculate y = x2 exp(x), we could put the following:
in an m-ile f12 .m (the function name and the m-ile name need to be the same). In the irst 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.,
|
|
|
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: |
|
num _points = 1000; h = logspace ( - 20 , - 1 , num _points ); for 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 ’);
|
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-ile f13 .m, then you need to use the @ operator to produce a function handle for quad:
|
|
|
If you used an inline function for f15, it’s already a function handle, and
|
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 use quad instead of integral; quad will return a 2nd output, which is the number of points it used.
