# Matlab functions

All these functions are released under the BSD-new license.

If you find any errors or can suggest some optimisations, let me know!

Note: I will try to keep the function headers on this page in sync with the Matlab files but I can't guarantee it.

## Miscellaneous

misc.zip: all the miscellaneous m-files.

expanddim.m

  function MM = expanddim(M, dim)

Take an interpolation or differentiation matrix and expand it so that it
works on multi-dimensional data sets.

Written by David A.W. Barton (david.barton@bristol.ac.uk) 2008

jacobian.m

  function J = jacobian(fun, sol, [h])

Construct a central difference approximation to the Jacobian of function

Default h = 1e-6

Written by David A.W. Barton (david.barton@bristol.ac.uk) 2008

simplespectrogram.m: requires windowing routines.

  function S = simplespectrogram(data, window, overlap, shift)

Calculate a short-time FFT of the data using overlapping Hamming windows.
Data is truncated if it is not correctly related to the window size.

data     Data to process (should be a 1d vector)
window   Width of the windows to use (default length(data)/8)
overlap  Size of the overlap between windows (default 50%)
shift    Whether to apply fftshift to the FFT results (default true)

For best performance, the window size should be a power of two. Result is
a matrix where the columns are the frequency direction and the rows are
the time direction.

Written by David A.W. Barton (david.barton@bristol.ac.uk) 2009

## Numerical integration

integration.zip: all the numerical integration m-files. Some of these require the polynomial and miscellaneous routines as well.

odegauss.m

  function sol = odegauss(ncoll, func, trange, y0, odeopt)

Integrate an ODE with an implicit Runge-Kutta collocation method at Gauss
points (order 2, 4 or 6). The (potentially nonlinear) system is solved
using a damped Newton iteration.

ncoll     Number of collocation points (1, 2 or 3)
func      Right-hand-side of the ODE in the form func(t, y)
trange    Start and end times for integration. (If integration output
at specific time points is required, use ODEGAUSSEVAL
afterwards.)
y0        Starting point for integration.
odeopt    Options structure returned by ODESET. Note: only the
InitialStep, Jacobian and Mass properties are used.

Written by David A.W. Barton (david.barton@bristol.ac.uk) 2011

odegausseval.m

  function y = odegausseval(sol, x)

Evaluate the solution obtained with ODEGAUSS at new time points.

sol      Previously obtained solution.
x        New time points.

Written by David A.W. Barton (david.barton@bristol.ac.uk) 2011

odegauss2.m: calls ODEGAUSS with ncoll = 1 (i.e., 2nd order implicit mid-point method)

odegauss4.m: calls ODEGAUSS with ncoll = 2 (i.e., 4th order)

odegauss6.m: calls ODEGAUSS with ncoll = 3 (i.e., 6th order)

## Polynomials

polynomials.zip: all the polynomial m-files.

bary_diff_matrix.m

  function [D M] = bary_diff_matrix(xnew, xbase, w)

Calculate the both the derivative and plain Lagrange interpolation matrix
using Barycentric formulae from Berrut and Trefethen, SIAM Review, 2004.

xnew     Interpolation points
xbase    Base points for interpolation
w        Weights calculated for base points

Written by David A.W. Barton (david.barton@bristol.ac.uk) 2008, 2011

bary_interp_matrix.m

  function M = bary_interp_matrix(xnew, xbase, w)

Calculate the Lagrange interpolation matrix using Barycentric formulae
from Berrut and Trefethen, SIAM Review, 2004.

xnew     Interpolation points
xbase    Base points for interpolation
w        Weights calculated for base points

Written by David A.W. Barton (david.barton@bristol.ac.uk) 2008, 2011

bary_weights_cheb1.m

  function w = bary_weights_cheb1(n)

Calculate Barycentric weights for the n degree Chebyshev polynomial of
the first kind.

Uses explicit formulae. See Berrut and Trefethen, SIAM Review, 2004.

Written by David A.W. Barton (david.barton@bristol.ac.uk) 2009

bary_weights_cheb2.m

  function w = bary_weights_cheb2(n)

Calculate Barycentric weights for the n degree Chebyshev polynomial of
the second kind.

Uses explicit formulae. See Berrut and Trefethen, SIAM Review, 2004.

Written by David A.W. Barton (david.barton@bristol.ac.uk) 2009

bary_weights_equispaced.m

  function w = bary_weights_equispaced(n)

Calculate Barycentric weights for the n degree polynomial with base
points equispaced.

Uses explicit formulae. See Berrut and Trefethen, SIAM Review, 2004.

Written by David A.W. Barton (david.barton@bristol.ac.uk) 2009

bary_weights.m

  function w = bary_weights(x)

Calculate Barycentric weights for the base points x.

Note: this algorithm may be numerically unstable for high degree
polynomials (e.g. 15+). If using linear or Chebyshev spacing of base
points then explicit formulae should then be used. See Berrut and
Trefethen, SIAM Review, 2004.

Written by David Barton (david.barton@bristol.ac.uk) 2008

poly_legendre.m

  function P = poly_legendre(n)

Return the n-th Legendre polynomial.

Written by David A.W. Barton (david.barton@bristol.ac.uk) 2009, 2010

poly_legendre_roots.m

  function root = poly_legendre_roots(n)

Return the roots of the n-th degree Gauss-Legendre polynomial on the
interval [-1, 1].

Written David A.W. Barton (david.barton@bristol.ac.uk) 2009, 2010

  function [w, x] = quad_gauss_legendre(n)

Quadrature weights for integrating f(x) over [-1, +1] at the Gauss nodes
(i.e., the roots of the appropriate Legendre polynomial).

Written by David A.W. Barton (david.barton@bristol.ac.uk) 2010, 2011

## Windowing

windowing.zip: all the windowing m-files.

hammingwindow.m

  function w = hammingwindow(n)

Return a Hamming window of width n as given by
w(i) = 0.54 - 0.46*cos(2*pi*i/(n - 1))
for 0 <= i <= n - 1.

Written by David A.W. Barton (david.barton@bristol.ac.uk) 2009

hannwindow.m

  function w = hannwindow(n)

Return a Hann window of width n as given by
w(i) = 0.5 - 0.5*cos(2*pi*i/(n - 1))
for 0 <= i <= n - 1.

Written by David A.W. Barton (david.barton@bristol.ac.uk) 2009