# 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 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)$$. numpy.ndarray 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(). tuple 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. numpy.ndarray 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. numpy.ndarray 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. numpy.ndarray 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). numpy.ndarray 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(). numpy.ndarray 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. numpy.ndarray The updated solution value.