RKtools

Some general utilities for analyzing Runge-Kutta methods.

Some of the routines expect as input a structured array \(rk\). This structure must have the fields \(A, b, c\), containing its Butcher coefficients. Optionally, it may represent an additive Runge-Kutta method or an embedded pair in which case it should also have \(Ahat\), \(bhat\), \(chat\) containing the coefficients of the secondary method.

internal_stab_explicit_butcher

function [stability] = internal_stab_explicit_butcher(A,b,c,spectrum,one_step_dt,p)

This function computes and plots both intermediate and one-step internal stability vector of an explicit Runge-Kutta scheme given its Butcher tableau.

Note that for an explicit Runge-Kutta scheme the stability functions are polynomials in the complex variable z.

Construct the intermediate stability functions psi_j (where j is the index of the stage).

Note that for an explicit scheme the intermediate stability polynomial associated to the first stage is always 1, i.e. psi_1 = 1. Therefore we just compute and plot the remaining (s-1) intermediate stability polynomials plus the one-step stability polynomial of the Runge-Kuatta method.

plotstabreg_func

function [contour_matrix] = plotstabreg_func(p,q,bounds,ls,lw)

plot the absolute stability region of a one-step method, given the stability function

Inputs:
  • p: coefficients of the numerator of the stability function

  • q: coefficients of the denominator of the stability function

if q is omitted, it is assumed that the function is a polynomial

Remaining inputs are optional:
  • bounds: bounds for region to compute and plot (default [-9 1 -5 5])

  • ls: line style (default ‘-r’)

  • lw: line width (default 2)

plotstabreg

function [contour_matrix] = plotstabreg(rk,plotbounds,ls,lw)

Plots the absolute stability region of a Runge-Kutta method, given the Butcher array

optimal_shuosher_form

function [v,alpha,beta] = optimal_shuosher_form(A,b,c)

L2_timestep_poly

function c = L2_timestep_poly(sdisc,p,q,eps,tol)

Find the absolutely timestep for a given combination of linear spatial discretization and stability function.

Also (optionally) plots the region of absolute stability and the eigenvalues.

The timestep is determined to within accuracy eps (default 10^-4).

The spectral stability condition is checked to within tol (default 10^-13).

semispectrum

function L = semispectrum(method,order,doplot,nx,cfl)

Plot spectra of various semi-discretizations of the advection equation

Current choices for method:
  • ‘fourier’: Fourier spectral method

  • ‘chebyshev’: Chebyshev spectral method

  • ‘updiff’: Upwind difference operators (linearized WENO)

  • ‘DG’: Discontinuous Galerkin method

The value of order matters only for the ‘updiff’ and ‘DG’ methods and selects the order of accuracy in those cases.

rk_stabfun

function [p,q] = rk_stabfun(rk)

Outputs the stability function of a RK method. The Butcher coefficients should be stored in rk.A, rk.b.

p contains the coefficients of the numerator

q contains the coefficients of the denominator

am_radius

function r = am_radius(A,b,c,eps,rmax)

Evaluates the Radius of absolute monotonicity of a Runge-Kutta method, given the Butcher array.

For an \(m\)-stage method, \(A\) should be an \(m \times m\) matrix and \(b\) and \(c\) should be column vectors of length m.

Accuracy can be changed by modifying the value of eps (default \(10^{-10}\)) Methods with very large radii of a.m. (>50) will require the default value of rmax to be increased.

The radius of absolute monotonicity is the largest value of \(r\) such that

\begin{align*} K(I+rA)^{-1} \ge & 0 \\ rK(I+rA)^{-1}e_m \le & e_{m+1} \end{align*}

where $$ K = \left(\begin{array}{c} A \\ b^T \end{array}\right) $$