Computing and Numerical Methods 1 Part 2: Numerical Methods Coursework: Bessel functions
Computing and Numerical Methods 1
Part 2: Numerical Methods
Completing this live file
You will submit only this self-contained completed livefile.No other scripts or other files will be marked. When completing this livefile:
· Please make sure all of your code is formatted as code using the "code" button on the Live Editor ribbon;
· Please make sure anything that is not code or code comments is formatted as text;
· Marks for each section and sub-section are shown in bold square brackets;
· Ensure that the entire livefile runs without errors or warnings before submitting!
Code
All code you write must have a comment clearly explaining what each line does. Comments should be brief but clear.
All local functions must:
· be entered in the Functions section at the bottom of this livefile;
· be named correctly according to the guidance below;
· have a comment on each line of code explaining its function and how it relates to the equations in the handouts or other theory;
· NOT call built-in MATLAB functions from outside the standard MATLAB toolbox C:\Program Files\MATLAB\R2022b\toolbox\matlab (but you may use these to check your functions).
To check the toolbox of a MATLAB function use the which command. The short name of the toolbox is shown after "toolbox\" For example:
which spline
returns C:\Program Files\MATLAB\R2022b\toolbox\matlab\polyfun\spline.m which means it is in the standard MATLAB toolbox, whereas
which fit
returns C:\Program Files\MATLAB\R2022b\toolbox\curvefit\curvefit\fit.m which means it is in the curvefit toolbox, and not the standard MATLAB toolbox.
Plots
You should ensure any plots you are asked to produce have:
· axis labels;
· x and y grid;
· all lines/markers clearly differentiable through use of colour or line/marker type;
· legend (where more than one dataset is plotted), placed so as to minimise overlap with data lines/markes;
· any spaced discrete point data (such as that in the table below) plotted as points (not lines);
· when plotting overlapping functions use solid, dashed, chain-dashed and dotted linetypes to maximise visibility of all lines.
Note that any plots that are not formatted as above (and exemplified in the Computer Labs Code livefiles) will not be marked.
Numerical outputs
Anywhere you are asked to output values you should:
· Use one of the appropriate built-in MATLAB functions e.g. sprintf (as in the Computer Lab Codes);
· Include a clear description of what the output variable represents;
· Absolute errors should be presented to 3 significant figures.
Submission
· You should submit your livefile on Blackboard by Friday 24th March 2023 at 1700.
· Please rename the file in the format NM1_CW_XXXXXX.mlx where XXXXXX is your CID number.
· As usual, please ensure you receive a confirmation email following submission. Only confirmed submitted files will be marked (partially complete Blackboard submissions will be ignored).
Introduction
Bessel functions describe electric fields, mechanical vibrations, heat conduction, optical diffraction and other phenomena involving cylindrical or spherical symmetry, for example the displacement of vibrating circular membranes, such as pressure transducer sensing elements and drum skins, as well as temperature distributions in heated cylindrical objects, such as combustor temperature probes and hot dogs in boiling water. The image below (source: https://commons.wikimedia.org/wiki/File:Vibrating_drum_Bessel_function.gif ) shows an animation of a vibrating drum skin. Everything you could possibly want to know about Bessel functions can be found at the DLMF: NIST Digital Library of Mathematical Functions, which is based on a famous mathemtical handbook known as "Abramowitz and Stegun", available as a pdf here.
There are various different kinds of Bessel functions, including the first kind , the second kind and the third kind . For each kind there are also different orders, e.g. integers The Bessel functions of the first kind for the first five integer orders (n = 0 to 4) are plotted below, using the built-in MATLAB function besselj. They look a bit like decaying sinusoids, except for some subtle differences as we shall see.
To see how this plot was made run the folllowing code:
openExample('matlab/PlotBesselFunctionsOfFirstKindExample')
Interestingly, the vertical displacement of the drumskin above is given by in polar coordinates. You should be able to see that, at a certain moment in time (e.g. at the limits of the motion), the vertical displacement plotted in the radial direction resembles 
in the plot above. You will learn how to obtain the full analytic solution to this sort of time-dependent differential equation later on in the undergraduate course.
In this coursework you are going to investigate Bessel functions of the first kind, , over the interval 
.
Part 1: Interpolation [9 marks]
Bessel functions are usually not amenable to straightforward evaluation (e.g. by hand or with a standard scientific calculator) and are often compiled in standard mathematical tables. One such table is shown below:
BesselTable= array2table([...
0 0
0.5 0.24227
1 0.44005
1.5 0.55794
2 0.57672
2.5 0.49709
3 0.33906
3.5 0.13738
4 -0.066043
4.5 -0.23106
5 -0.32758
5.5 -0.34144
6 -0.27668
6.5 -0.15384
7 -0.0046828
7.5 0.13525
8 0.23464
8.5 0.27312
9 0.24531
9.5 0.16126
10 0.043473...
], "VariableNames",["x","J"])
BesselTable = 21×2 table
| x | J |
1 | 0 | 0 |
2 | 0.5000 | 0.2423 |
3 | 1 | 0.4400 |
4 | 1.5000 | 0.5579 |
5 | 2 | 0.5767 |
6 | 2.5000 | 0.4971 |
7 | 3 | 0.3391 |
8 | 3.5000 | 0.1374 |
9 | 4 | -0.0660 |
10 | 4.5000 | -0.2311 |
11 | 5 | -0.3276 |
12 | 5.5000 | -0.3414 |
13 | 6 | -0.2767 |
14 | 6.5000 | -0.1538 |
Investigate the different methods of interpolation covered in the course using the interpolant data in the above table as follows:
· M1.1: Write a local MATLAB function called IntPoly that fits a polynomial of the maximum possible order to the tabular data using the first approach in Handout 4 involving inversion of a matrix of x-values. [1]
· M1.2: Write a local MATLAB function called LagrangePoly to compute the interpolating Lagrange polynomial that passes through all inerpolant data points. HINT: You may find the matlab functions poly and polyval useful. [2]
· M1.3: Write a local MATLAB function called IntSpline that performs cubic spline interpolation with natural (free-runout) end conditions. [1]
All functions should:
· take as input a vector of query points that is evenly spread across the interval and spaced by 0.1
· output the evaluation of the interpolation function at the query points.
Evaluate the performance of your local functions as follows: [1]
· E1.1: Plot the interpolated data (as a continuous curve) alongside the tabular data (as points) for all methods on a single plot.
· E1.2: Evaluate and plot the absolute error in the interpolation compared to besselj. Plot the absolute error as a function of x for all methods on a single plot.
· E1.3: Compute and output the time required to perform the interpolation along with the 
-norm of the absolute error across the interval using e.g. sprintf (as in the Computer Lab Codes) to clearly state which values relate to which method.
Questions
· Q1.1: How does the performance of the various methods compare in terms of accuracy and computation time? [1]
· Q1.2: Methods 1. and 2. should generate equivalent polynomials according to Handout 4. Do they? Try to explain any discrepancies. [1]
· Q1.3: What happens to the overall error (-norm) in method 2 if, instead of using the tabular data, you interpolate data sampled directly from besselj with the same number of points (i.e. 21) but using Chebyshev spacing? [1]
· Q1.4: Compare method 3. with the built-in MATLAB function spline. Are the outputs the same? Try to explain any discrepancies. [1]
Hints
Method M1.2: One way to approach this is to use a loop to step through the 
and for each one create a matrix of basis polynomial coefficients, using the built-in MATLAB function poly to find the coefficients of the (basis) polynomial whose roots are all the except one. You can then multiply this by a vector of 
in order to obtain the scaled basis polynomials, and then use built-in MATLAB function polyval to evaluate the resulting polynomial coefficients for the query points.
Part 1: Your Solution
% Enter your code here and your local functions at the bottom.
Enter explanation here.
Part 2: Numerical Integration [8 marks]
For integer values of n, the bessel function of the first kind may be defined in integral form using Bessel's integrals as:
Investigate the different methods of numerical integration covered in the course to evaluate Bessel's integral for the first order Bessel Function of the first kind, as follows:
· M2.1: Write a local MATLAB function called trapzJ1 that uses the Trapezium rule. [1]
· M2.2: Write a local MATLAB function called simpsonJ1 that uses Simpson's 1/3 rule. [1]
· M2.3: Write a local MATLAB function called gaussJ1 that uses Gauss-Legendre quadrature. [3]
All functions should:
· take as input a vector x that is evenly spread across the interval and spaced by 0.1, as well as the number of points used in the integration, N.
· output 
.
You should evaluate the performance of your local functions as follows [2]:
· E2.1: Plot the approximated function alongside values computed using besselj for all methods on a single plot.
· E2.2: Evaluate and plot the absolute error compared to besselj as a function of x for all methods on a single plot.
· E2.3: On a single figure, plot the -norm of the absolute error across the interval vs the number of points, N, used in the integration for each method for a sensible range of N. Plot a horizontal line representing double-precision accuracy, and for each method work out how many points are required to evaluate to within machine (double) precision of besselj.
· E2.4: Write a simple loop to evaluate the integral for the vector x using the built-in MATLAB function integral. Plot the error across the interval and evaluate the computation time required to achieve double-precision accuracy.
Questions
· Q2.1: How does the computation time of the three methods M2.1,M2.2, M2.3 and the looped integral function compare? Do you think the function besselj uses numerical integration to evaluate Bessel functions? [1]
Hints
Method M2.1: Use the built-in MATLAB function trapz.
Method M2.2: Use an anonymous function (as in the Computer Lab Codes), that takes x and tau as the inputs, in order to define the integrand.
Method M2.3: You will need to evaluate the Legendre polynomials to higher orders than in the Class Handouts.
One way to approach this is to use the MATLAB functions conv and polyder to evaluate the definition of
in Handout 3:
The 
term can be evaluated using the built-in MATLAB function conv, which multiplies two polynomials and returns the coefficients of the resulting polynomial (a process known as convolution):
pn = [1 0 -1]; % Polynomial coefficients for polynomial x^2-1
res = 1; % Initialise the result
N = 3; % Power to which you want to calculate the input polynomial
for k1 = 1:N
res = conv(res, pn); % conv . The first iteration will simply multiply 1 and [1 0 1] to obtain [1 0 1].
end
% res will now be the coefficients of (x^2-1)^N expanded as a polynomial.
The coefficients of the derivative of a polynomial can be found using built-in MATLAB function polyder. You will need to code a loop similar to the one for conv to calculate the Nth derivative. Yo can also use polyder and polyval to evaluate the weights.
To obtain the roots of the polynomial based on the coefficient simply use the built-in MATLAB function roots. To evaluate a polynomial for a given vector of input values simply use the built-in MATLAB function polyval.
Another way to generate the Legendre polynomials is using Bonnet's recursion formula:
Starting with and this allows all subsequent polynomials to be computed. For example, is found from:
.
Based on this it is possible to code a recurrence relation to calculate the ceofficients of any Legendre polynomial. To illustrate this consider the coefficients of the first few polynomials:
L0_coeffs = 1
L1_coeffs = [1 0] = [L0_coeffs 0]/1
L2_coeffs = 1/2*[3 0 -1] = ( 3*[L1_coeffs 0] - [0 0 L0coeffs] )/2
You should be able to spot a pattern forming that corresponds to Bonnet's recursion formula. Obviously you will need to use a loop and work out how to increment the number of coefficients for each iteration. You can check your polynomials against the symbolic math toolbox function legendreP.
Part 2: Your Solution
% Enter your code here and your local functions at the bottom.
Enter explanation here.
Part 3: Function Roots [8 marks]
Compare the various root finding methods that are amenable to this function, without knowing its derivative, to find all roots on the interval. Compare with sine and cosine. Compare with fzero.
Investigate the different methods of root finding covered in the course to find all of the roots on the interval for the zeroth order Bessel function of the first kind, , as follows:
· M3.1: Write a local MATLAB function called fzeroJ that uses the built-in MATLAB function fzero. [1]
· M3.2: Write a local MATLAB function called bisectionJ that uses the bisection method. [1]
· M3.3: Write a local MATLAB function called newtonJ that uses Newton's method. [1]
· M3.4: Write a local MATLAB function called secantJ that uses the Secant method. [1]
All functions should:
· take as input the order of the Bessel function and the start and end points of the interval (i.e. [0 10]).
· output the roots to a tolerance of 1e-6 (in other words stop iterating when the solution is changing by less than the tolerance), if this is possible within a computation time of say 10 minutes (use Ctrl+C or Command+. to stop the computation if not).
· NOT use predetermined search intervals, but rather determine these from the function itself (so that in principle your code would function for any order/kind of Bessel function and/or outside the interval).
· all use the same method for determining the search interval.
You should evaluate the performance of your local functions as follows: [1]
· E3.1: Plot the approximated roots alongside the function evaluated using besselj, using a different plot for each method.
· E3.2: Compute and output the time required to compute all roots on the interval using e.g. sprintf (as in the Computer Lab Codes) to clearly state which values relate to which method.
Questions
· Q3.1: How do the various methods compare in terms of accuracy and computation time? [1]
· Q3.2: Are your roots to within machine precision of the true values on page 409 of Abramowitz and Stegun? If not, why not? [1]
· Q3.3: Does your function work for other orders of the Bessel function of the first kind, , and/or for roots outside the interval (i.e. )? If so, when does it stop working and why? How would you need to modify your code to get it to work? N.B. you should not try to modify your function to get it to work on other functions or outside the interval, as this problem is hard enough already, but you should only try to understand (and explain) why it doesn't work. [1]
Hints
One way to find the search interval is to start from the origin and find the first data point where the sign has changed, then search between that point and the next. You can then continue this to the next root.
Method 3.3: You will need to work out the derivative. This can be found from the series definition of the function, and you can also look it up at DLMF: §10.6 Recurrence Relations and Derivatives ‣ Bessel and Hankel Functions ‣ Chapter 10 Bessel Functions (nist.gov).
Part 3: Your Solution
% Enter your code here and your local functions at the bottom.
Enter explanation here.
Your Functions
Part 1:
% Local functions for part 1.
Part 2:
% Local functions for part 2.
Part 3:
% Local functions for part 3.
