assignment1.dg1 module

Module for solving a 1D conservation law via DG.

Adapted from a Discontinuous Galerkin (DG) solver written by Per Olof-Persson.

Check out an example notebook using these utilities to solve the problem.

class assignment1.dg1.DG1Solver(num_intervals, p_order, total_time, dt, get_initial_data=None, points_on_ref_int=None)[source]

Bases: object

Discontinuous Galerkin (DG) solver for the 1D conservation law.

\[\frac{\partial u}{\partial t} + \frac{\partial u}{\partial x} = 0\]

on the interval \(\left[0, 1\right]\). By default, uses Gaussian-like initial data

\[u(x, 0) = \exp\left(-\left(\frac{x - \frac{1}{2}}{0.1} \right)^2\right)\]

but \(u(x, 0)\) can be specified via get_initial_data.

We represent our solution via the \((p + 1) \times n\) rectangular matrix:

\[\begin{split}\mathbf{u} = \left[ \begin{array}{c c c c} u_0^1 & u_0^2 & \cdots & u_0^n \\ u_1^1 & u_1^2 & \cdots & u_1^n \\ \vdots & \vdots & \ddots & \vdots \\ u_p^1 & u_p^2 & \cdots & u_p^n \end{array}\right]\end{split}\]

where each column represents one of \(n\) sub-intervals and each row represents one of the \(p + 1\) node points within each sub-interval.

Parameters:
  • num_intervals (int) – The number \(n\) of intervals to divide \(\left[0, 1\right]\) into.
  • p_order (int) – The degree of precision for the method.
  • total_time (float) – The amount of time to run the solver for (starts at \(t = 0\)).
  • dt (float) – The timestep to use in the solver.
  • get_initial_data (callable) – (Optional) The function to use to evaluate \(u(x, 0)\) at the points in our solution. Defaults to get_gaussian_like_initial_data().
  • points_on_ref_int (function) – (Optional) The method used to partition the reference interval \(\left[0, h\right]\) into node points. Defaults to get_evenly_spaced_points().
get_mass_and_stiffness_matrices()[source]

Get the mass and stiffness matrices for the current solver.

Uses pre-computed mass matrix and stiffness matrix for \(p = 1\), \(p = 2\) and \(p = 3\) degree polynomials and computes the matrices on the fly for larger \(p\).

Depends on the sub-interval width h and the order of accuracy p_order.

Return type:tuple
Returns:Pair of mass and stiffness matrices, both with p_order + 1 rows and columns.
ode_func(u_val)[source]

Compute the right-hand side for the ODE.

When we write

\[\begin{split}M \dot{\mathbf{u}} = K \mathbf{u} + \left[ \begin{array}{c c c c c} u_p^2 & u_p^3 & \cdots & u_p^n & u_p^1 \\ 0 & 0 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & 0 & 0 \\ -u_p^1 & -u_p^2 & \cdots & -u_p^{n-1} & -u_p^n \\ \end{array}\right]\end{split}\]

we specify a RHS \(f(u)\) via solving the system.

Parameters:u_val (numpy.ndarray) – The input to \(f(u)\).
Return type:numpy.ndarray
Returns:The value of the slope function evaluated at u_val.
update()[source]

Update the solution for a single time step.

We use ode_func() to compute \(\dot{u} = f(u)\) and pair it with an RK method (low_storage_rk()) to compute the updated value.

class assignment1.dg1.MathProvider[source]

Bases: object

Mutable settings for doing math.

For callers that wish to swap out the default behavior – for example, to use high-precision values instead of floats – this class can just be monkey patched on the module.

The module assumes through-out that solution data is in NumPy arrays, but the data inside those arrays may be any type.

Note

The callers assume exp_func is a vectorized exponential that can act on a NumPy array containing elements of the relevant type.

Note

The zeros constructor should also be able to take the order argument (and should produce a NumPy array).

exp_func
linspace
mat_inv
num_type

alias of __builtin__.float

solve
zeros
assignment1.dg1.find_matrices(p_order, points_on_ref_int=None)[source]

Find mass and stiffness matrices.

We do this on the reference interval \(\left[-1, 1\right]\). By default we use the evenly spaced points

\[x_0 = -1, x_1 = -(p - 2)/p, \ldots, x_p = 1\]

but the set of nodes to use on the reference interval can be specified via the points_on_ref_int argument. With our points, we compute the polynomials \(\varphi_j(x)\) such that \(\varphi_j\left(x_i\right) = \delta_{ij}\). We do this by writing

\[\varphi_j(x) = \sum_{n = 0}^p c_n^{(j)} L_n(x)\]

where \(L_n(x)\) is the Legendre polynomial of degree \(n\). With this representation, we need to solve

\[\begin{split}\left[ \begin{array}{c c c c} L_0(x_0) & L_1(x_0) & \cdots & L_p(x_0) \\ L_0(x_1) & L_1(x_1) & \cdots & L_p(x_p) \\ \vdots & \vdots & \ddots & \vdots \\ L_0(x_p) & L_1(x_p) & \cdots & L_p(x_p) \end{array}\right] \left[ \begin{array}{c c c c} c_0^{(0)} & c_0^{(1)} & \cdots & c_0^{(p)} \\ c_1^{(0)} & c_1^{(1)} & \cdots & c_1^{(p)} \\ \vdots & \vdots & \ddots & \vdots \\ c_p^{(0)} & c_p^{(1)} & \cdots & c_p^{(p)} \end{array}\right] = \left(\delta_{ij}\right) = I_{p + 1}\end{split}\]

Then use these to compute the mass matrix

\[M_{ij} = \int_{-1}^1 \varphi_i(x) \varphi_j(x) \, dx\]

and the stiffness matrix

\[K_{ij} = \int_{-1}^1 \varphi_i'(x) \varphi_j(x) \, dx\]

Utilizing the fact that

\[\left\langle L_n, L_m \right\rangle = \int_{-1}^1 L_n(x) L_m(x) \, dx = \frac{2}{2n + 1} \delta_{nm}\]

we can compute

\[M_{ij} = \left\langle \varphi_i, \varphi_j \right\rangle = \sum_{n, m} \left\langle c_n^{(i)} L_n, c_m^{(j)} L_m \right\rangle = \sum_{n = 0}^p \frac{2}{2n + 1} c_n^{(i)} c_n^{(j)}.\]

Similarly

\[\begin{split}\left\langle L_n'(x), L_m(x) \right\rangle = \begin{cases} 2 & \text{ if } n > m \text{ and } n - m \equiv 1 \bmod{2} \\ 0 & \text{ otherwise}. \end{cases}\end{split}\]

gives

\[\begin{split}\begin{align*} K_{ij} &= \left\langle \varphi_i', \varphi_j \right\rangle = \sum_{n, m} \left\langle c_n^{(i)} L_n', c_m^{(j)} L_m \right\rangle \\ &= 2 \left(c_0^{(j)} \left(c_1^{(i)} + c_3^{(i)} + \cdots\right) + c_1^{(j)} \left(c_2^{(i)} + c_4^{(i)} + \cdots\right) + \cdots + c_{p - 1}^{(j)} c_p^{(i)}\right) \\ \end{align*}\end{split}\]

(For more general integrals, one might use Gaussian quadrature. The largest degree integrand \(\varphi_i \varphi_j\) has degree \(2 p\) so this would require \(n = p + 1\) points to be exact up to degree \(2(p + 1) - 1 = 2p + 1\).)

Parameters:
  • p_order (int) – The degree of precision for the method.
  • points_on_ref_int (function) – (Optional) The method used to partition the reference interval \(\left[0, h\right]\) into node points. Defaults to get_evenly_spaced_points().
Return type:

tuple

Returns:

Pair of mass and stiffness matrices, square numpy.ndarray of dimension p_order + 1.

assignment1.dg1.gauss_lobatto_points(start, stop, num_points)[source]

Get the node points for Gauss-Lobatto quadrature.

Using \(n\) points, this quadrature is accurate to degree \(2n - 3\). The node points are \(x_1 = -1\), \(x_n = 1\) and the interior are \(n - 2\) roots of \(P'_{n - 1}(x)\).

Though we don’t compute them here, the weights are \(w_1 = w_n = \frac{2}{n(n - 1)}\) and for the interior points

\[w_j = \frac{2}{n(n - 1) \left[P_{n - 1}\left(x_j\right)\right]^2}\]

This is in contrast to the scheme used in Gaussian quadrature, which use roots of \(P_n(x)\) as nodes and use the weights

\[w_j = \frac{2}{\left(1 - x_j\right)^2 \left[P'_n\left(x_j\right)\right]^2}\]

Note

This method is not generic enough to accommodate non-NumPy types as it relies on the numpy.polynomial.legendre.

Parameters:
  • start (float) – The beginning of the interval.
  • stop (float) – The end of the interval.
  • num_points (int) – The number of points to use.
Return type:

numpy.ndarray

Returns:

1D array, the interior quadrature nodes.

assignment1.dg1.get_evenly_spaced_points(start, stop, num_points)[source]

Get points on an interval that are evenly spaced.

This is intended to be used to give points on a reference interval when using DG on the 1D problem.

Parameters:
  • start (float) – The beginning of the interval.
  • stop (float) – The end of the interval.
  • num_points (int) – The number of points to use on the interval.
Return type:

numpy.ndarray

Returns:

The evenly spaced points on the interval.

assignment1.dg1.get_gaussian_like_initial_data(node_points)[source]

Get the default initial solution data.

In this case it is

\[u(x, 0) = \exp\left(-\left(\frac{x - \frac{1}{2}}{0.1} \right)^2\right)\]
Parameters:node_points (numpy.ndarray) – Points at which evaluate the initial data function.
Return type:numpy.ndarray
Returns:The \(u\)-values at each node point.
assignment1.dg1.get_legendre_matrix(points, max_degree=None)[source]

Evaluate Legendre polynomials at a set of points.

If our points are \(x_0, \ldots, x_p\), this computes

\[\begin{split}\left[ \begin{array}{c c c c} L_0(x_0) & L_1(x_0) & \cdots & L_d(x_0) \\ L_0(x_1) & L_1(x_1) & \cdots & L_d(x_p) \\ \vdots & \vdots & \ddots & \vdots \\ L_0(x_p) & L_1(x_p) & \cdots & L_d(x_p) \end{array}\right]\end{split}\]

by utilizing the recurrence

\[n L_n(x) = (2n - 1) x L_{n - 1}(x) - (n - 1) L_{n - 2}(x)\]
Parameters:
  • points (numpy.ndarray) – 1D array. The points at which to evaluate Legendre polynomials.
  • max_degree (int) – (Optional) The maximum degree of Legendre polynomial to use. Defaults to one less than the number of points (which will produce a square output).
Return type:

numpy.ndarray

Returns:

The 2D array containing the Legendre polynomials evaluated at our input points.

assignment1.dg1.get_node_points(num_points, p_order, step_size=None, points_on_ref_int=None)[source]

Return node points to split unit interval for DG.

Parameters:
  • num_points (int) – The number \(n\) of intervals to divide \(\left[0, 1\right]\) into.
  • p_order (int) – The degree of precision for the method.
  • step_size (float) – (Optional) The step size \(1 / n\).
  • points_on_ref_int (function) – (Optional) The method used to partition the reference interval \(\left[0, h\right]\) into node points. Defaults to get_evenly_spaced_points().
Return type:

numpy.ndarray

Returns:

The \(x\)-values for the node points, with p_order + 1 rows and \(n\) columns. The columns correspond to each sub-interval and the rows correspond to the node points within each sub-interval.

assignment1.dg1.low_storage_rk(ode_func, u_val, dt)[source]

Update an ODE solutuon with an order 2/4 Runge-Kutta function.

The method is given by the following Butcher array:

\[\begin{split}\begin{array}{c | c c c c} 0 & 0 & & & \\ 1/4 & 1/4 & 0 & & \\ 1/3 & 0 & 1/3 & 0 & \\ 1/2 & 0 & 0 & 1/2 & 0 \\ \hline & 0 & 0 & 0 & 1 \end{array}\end{split}\]

It is advantageous because the updates \(k_j\) can be over-written at each step, since they are never re-used.

One can see that this method is order 2 for general \(\dot{u} = f(u)\) by verifying that not all order 3 node conditions are satisfied. For example:

\[\frac{1}{3} \neq \sum_i b_i c_i^2 = 0 + 0 + 0 + 1 \cdot \left(\frac{1}{2}\right)^2\]

However, for linear ODEs, the method is order 4. To see this, note that the test problem \(\dot{u} = \lambda u\) gives the stability function

\[R\left(\lambda \Delta t\right) = R(z) = 1 + z + \frac{z^2}{2} + \frac{z^3}{6} + \frac{z^4}{24}\]

which matches the Taylor series for \(e^z\) to order 4.

See Problem Set 3 from Persson’s Math 228A for more details.

Parameters:
  • ode_func (callable) – The RHS in the ODE \(\dot{u} = f(u)\).
  • u_val (numpy.ndarray) – The input to \(f(u)\).
  • dt (float) – The timestep to use.
Return type:

numpy.ndarray

Returns:

The updated solution value.