COSC2500/COSC7500—Semester 2 2024R

Java Python COSC2500/COSC7500—Semester 2, 2024

Exercises—due 3:00 pm Friday 30th 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 For the polynomial f (x)  = x4 − 2x3, 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.4) and (1.6)).

a) Choose different step sizes hand 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.

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’twrite?

R1.2 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 size from that.

Note that 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’twrite?

R1.3 Integration over discontinuities. Integrate the function

from x = 0 to x = 3, for varying step sizes, using the trapezoidal method. Plot the error as a function of step size. A Matlab function calculating this function, f_r13a.m, is available on Blackboard.  Note that the integral is equal to 2e − 2.

Don’t just integrate from x = 1 to x = 2! We want to see what integrating over the discontinuities does!

Checklist

Did you integrate the function for varying step sizes? Did you explain how you did it?

Did you plot graphs showing the error?

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’twrite?

R1.4 Integration over an infinite interval. 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.

Checklist

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 fixed 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’twrite?

R1.5 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:

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 approxima- tion. You can use either equation (1.15) or trial-and-error and your own judgment to find the best order of polynomial to use.

Checklist

Did you fit a suitable function or functions to the data?

Did you include a graph comparing your fit to the data?

Did you explain how you chose the best fit?

Did you include a graph or fit measure values illustrating your choice?

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 6 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 Calculate the 2nd derivative and higher derivatives off (x) = 1 + x + x2 + x3 + x4 + x5 + x6 + x7 + x8 using the finite difference formulae for higher derivatives. How do the errors vary with step size?

A1.2 How would you calculate the integral

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 β .

A1.3 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.)  Compare your esti- mate of the error with actual error determined from comparison with the analytical solution, if an analytical solution is available.

A1.4 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.5 Compare polynomials and real Fourier series (i.e., sine–cosine series) for fitting the data in a15data . m. Matlab code for calculating/fitting a Fourier series is available on Blackboard:  sine_cosine_fit_int .m, which uses a Fourier integral, and sine_cosine_fit_ls.m, which uses a least-squares fit.

A1.6 Write a program to calculate the increase in the viscous drag coefficient due to a wall using the approximation given by Chaoui and Feuillebois (2003).

A1.7 How are the special functions available on a scientific calculator calculated? You might like to write some code to calculate these to test the efficiency and accuracy of the methods.

A1.8 Consider the algorithm for the calculation of Wigner 3jfunctions described by Brock (2001). Note that there are a number of special cases. How impor- tant is it to treat these cases separately? AMatlab implementation of Brock’s algorithm, wigner3j . mis 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.)

Programming hints

R1.1 Since you want to look at the effect of step size on error, you want to find 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 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.

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 = 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.2 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 use quad instead of integral; quad will return a 2nd output, which is the number of points it used.

R1.5 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