Spectral Element Method

\begin{align*} \begin{cases} u_t = a(t, x)u_{xx} + b(t, x) u_x + c(t, x) u \qquad \forall x\in [L, R], t\in[0, T]\\ u(t, L) = u_L(t)\\ u(t, R) = u_R(t)\\ u(0, x) = u_0(x) \end{cases}. \end{align*}

\begin{align*} \int_L^R u_t(t, x)\phi(x)\,dx =& \int_L^R a(t, x)u_{xx}(t, x)\phi(x)\,dx \nonumber\\ &+ \int_L^R b(t, x) u_x(t, x)\phi(x)\,dx \nonumber\\ &+ \int_L^R c(t, x) u(t, x)\phi(x)\,dx.\qquad\label{eq:weakformLR} \end{align*}

\begin{align*} \int_L^R f(t, x)\,dx = \frac{d}2 \int_{-1}^1 f(t, x(\xi))\,d\xi. \end{align*}

\begin{align*} \frac{d}2\int_{-1}^1 u_t(t, \xi)\phi(\xi)\,d\xi =& \frac2d\int_{-1}^1 a(t, \xi)u_{\xi\xi}(t, \xi)\phi(\xi)\,d\xi \nonumber\\ &+ \int_{-1}^1 b(t, \xi) u_\xi(t, \xi)\phi(\xi)\,d\xi \nonumber\\ &+ \frac{d}2\int_{-1}^1 c(t, \xi) u(t, \xi)\phi(\xi)\,d\xi.\qquad\label{eq:weakform} \end{align*}

\begin{align*} \int_{-1}^1 a u_{\xi\xi} \phi\,d\xi &= u_{\xi\xi}a\phi\Big|_{-1}^1 - \int_{-1}^1 u_{\xi}(a\phi)_{\xi}\,d\xi\\ &= u_{\xi\xi}a\phi\Big|_{-1}^1 - \int_{-1}^1 u_{\xi}a\phi^\prime\,d\xi - \int_{-1}^1 u_{\xi}a_{\xi}\phi\,d\xi\\ &= - \int_{-1}^1 u_{\xi}a\phi^\prime\,d\xi - \int_{-1}^1 u_{\xi}a_{\xi}\phi\,d\xi, \end{align*}

\begin{align*} \int_{-1}^1 a u_{\xi\xi} \phi\,d\xi = - \int_{-1}^1 u_{\xi}a\phi^\prime\,d\xi - \frac{d}2\int_{-1}^1 u_{\xi}a_x\phi\,d\xi. \end{align*}

\begin{align*} \frac{d}2\int_{-1}^1 u_t\phi\,d\xi = -\frac2d\int_{-1}^1 au_{\xi}\phi^\prime\,d\xi + \int_{-1}^1 \left(b-a_x\right) u_\xi\phi\,d\xi + \frac{d}2\int_{-1}^1 c u\phi\,d\xi.\label{eq:bypart} \end{align*}

\begin{align*} u(t, \xi) &\approx \sum_{p=0}^N U_p(t)\ell_p(\xi),\\ \phi(\xi) &\approx \sum_{i=0}^N \phi_i\ell_i(\xi),\\ \int_{-1}^1 f(\xi)\,d\xi &\approx \sum_{n=0}^N f(\xi_n)w_n, \end{align*}

\begin{align*} \frac{d}2~\sum_{p=0}^N\sum_{i=0}^N\sum_{n=0}^N \dot U_p(t)\ell_p(\xi_n)\phi_i\ell_i(\xi_n)w_n =& -\frac2d~ \sum_{p=0}^N\sum_{i=0}^N\sum_{n=0}^N U_p(t)\ell_p^\prime(\xi_n)a(t, \xi_n)\phi_i\ell_i^\prime(\xi_n)w_n \nonumber\\ &+~ \sum_{p=0}^N\sum_{i=0}^N\sum_{n=0}^N U_p(t)\ell_p^\prime(\xi_n)(b-a_x)(t, \xi_n)\phi_i\ell_i(\xi_n)w_n \nonumber\\ &+~ \frac{d}2~\sum_{p=0}^N\sum_{i=0}^N\sum_{n=0}^N U_p(t)\ell_p(\xi_n)c(t, \xi_n)\phi_i\ell_i(\xi_n)w_n \nonumber\\ \label{eq:fullapprox} \end{align*}

\begin{align*} \ell_p(\xi_n) = \begin{cases} 1 &\mbox{ if }p=n\\ 0 &\mbox{ if }p\ne n \end{cases}. \end{align*}

\begin{align*} \frac{d}2~\sum_{i=0}^N \dot U_i(t)\ell_i(\xi_i)\phi_i\ell_i(\xi_i)w_i = \frac{d}2~ \sum_{i=0}^N \dot U_i(t)\phi_iw_i. \end{align*}

\begin{align*} \frac{d}2~ \sum_{i=0}^N \dot U_i(t)\phi_iw_i =& -\frac2d~ \sum_{p=0}^N\sum_{i=0}^N\sum_{n=0}^N U_p(t)\ell_p^\prime(\xi_n)a(t, \xi_n)\phi_i\ell_i^\prime(\xi_n)w_n \nonumber\\ &+~ \sum_{p=0}^N\sum_{i=0}^N U_p(t)\ell_p^\prime(\xi_i)(b-a_x)(t, \xi_i)\phi_iw_i \nonumber\\ &+~ \frac{d}2~\sum_{i=0}^N U_i(t)c(t, \xi_n)\phi_iw_i \nonumber\\ =& -\frac2d~ \sum_{p=0}^N\sum_{i=0}^N\sum_{n=0}^N U_pD_{np}A_i\phi_iD_{ni}w_n \nonumber\\ &+~ \sum_{p=0}^N\sum_{i=0}^N U_pD_{ip}(b-a_x)_i\phi_iw_i \nonumber\\ &+~ \frac{d}2~\sum_{i=0}^N U_iC_i\phi_iw_i, \label{eq:phiieq} \end{align*}

\begin{align*} \frac{d}{2}~\dot U_i(t) w_i =& -\frac2d~ \sum_{p=0}^N\sum_{n=0}^N U_p(t)D_{np}A_iD_{ni}w_n \nonumber\\ &+~ \sum_{p=0}^N U_p(t)D_{ip}(b-a_x)_iw_i\nonumber\\ &+~ \frac{d}2~ U_i(t)C_iw_i\qquad \forall i=1, 2, \ldots, N-1. \label{eq:ODEsummation} \end{align*}

With the same notation as in (\ref{eq:ODEsummation}), we have

\(\sum_{p=0}^N\sum_{n=0}^N U_p(t)D_{np}A_i\phi_iD_{ni}w_n\) is the \(i\)-th element of the vector \(D^TAWD\vec U\);

\(\sum_{p=0}^N U_p(t)D_{ip}(b-a_x)_i\phi_iw_i\) is the \(i\)-th element of the vector \((B-A_x)WD\vec U\), where

:nbsphinx-math:`begin{align*} &vec U = left( begin{array}{c} U_0(t)\ U_1(t)\ U_2(t)\ vdots\ U_N(t) end{array} right), A_x = left( begin{array}{ccccc} a_x(t, xi_0) & & & & 0 \ & a_x(t, xi_1)& & & \ & & a_x(t, xi_2)& & \ & & & ddots & \ 0 & & & & a_x(t, xi_N) end{array} right), \\ &A = left( begin{array}{ccccc} a(t, xi_0) & & & & 0 \ & a(t, xi_1)& & & \ & & a(t, xi_2)& & \ & & & ddots & \ 0 & & & & a(t, xi_N) end{array} right), B = left( begin{array}{ccccc} b(t, xi_0) & & & & 0 \ & b(t, xi_1)& & & \ & & b(t, xi_2)& & \ & & & ddots & \ 0 & & & & b(t, xi_N) end{array} right), \\ &C = left( begin{array}{ccccc} c(t, xi_0) & & & & 0 \ & c(t, xi_1)& & & \ & & c(t, xi_2)& & \ & & & ddots & \ 0 & & & & c(t, xi_N) end{array} right), W = left( begin{array}{ccccc} w_0 & & & & 0 \ & w_1 & & & \ & & w_2 & & \ & & & ddots & \ 0 & & & & w_N end{array} right), \\ &D = left( begin{array}{ccccc}

D_{00} & D_{01} & D_{02} & cdots & D_{0N} \ D_{10} & D_{11} & D_{12} & cdots & D_{1N} \ D_{20} & D_{21} & D_{22} & cdots & D_{2N} \

vdots & vdots & vdots & ddots & vdots \

D_{N0} & D_{N1} & D_{N2} & cdots & D_{NN}

end{array} right) = left( begin{array}{ccccc}

ell_0^prime(xi_0) & ell_1^prime(xi_0) & ell_2^prime(xi_0) & cdots & ell_N^prime(xi_0) \ ell_0^prime(xi_1) & ell_1^prime(xi_1) & ell_2^prime(xi_1) & cdots & ell_N^prime(xi_1) \ ell_0^prime(xi_2) & ell_1^prime(xi_2) & ell_2^prime(xi_2) & cdots & ell_N^prime(xi_2) \

vdots & vdots & vdots & ddots & vdots \

ell_0^prime(xi_N) & ell_1^prime(xi_N) & ell_2^prime(xi_N) & cdots & ell_N^prime(xi_N)

end{array} right). end{align*}`

\begin{align*} \frac{d}2 W\dot{\vec{U}} = \left( -\frac2d D^TAWD + \left(B-A_x\right)WD + \frac{d}2CWD\right) \vec U. \end{align*}

\begin{align*} \left(\int_{x_0}^{x_1} + \int_{x_1}^{x_2} + \int_{x_2}^{x_3}\right) u_t(t, x)\phi(x)\,dx =& \left(\int_{x_0}^{x_1} + \int_{x_1}^{x_2} + \int_{x_2}^{x_3}\right) a(t, x)u_{xx}(t, x)\phi(x)\,dx \nonumber\\ &+ \left(\int_{x_0}^{x_1} + \int_{x_1}^{x_2} + \int_{x_2}^{x_3}\right) b(t, x) u_x(t, x)\phi(x)\,dx \nonumber\\ &+ \left(\int_{x_0}^{x_1} + \int_{x_1}^{x_2} + \int_{x_2}^{x_3}\right) c(t, x) u(t, x)\phi(x)\,dx.\qquad\qquad\label{eq:weakformLRelem} \end{align*}

\begin{align*} \int_{x_0}^{x_1} u_t(t, x)\phi^{(1)}(x)\,dx =& \int_{x_0}^{x_1} a(t, x)u_{xx}(t, x)\phi^{(1)}(x)\,dx \nonumber\\ &+ \int_{x_0}^{x_1} b(t, x) u_x(t, x)\phi^{(1)}(x)\,dx \nonumber\\ &+ \int_{x_0}^{x_1} c(t, x) u(t, x)\phi^{(1)}(x)\,dx,\label{eq:weakform1}\\ \int_{x_1}^{x_2} u_t(t, x)\phi^{(2)}(x)\,dx =& \int_{x_1}^{x_2} a(t, x)u_{xx}(t, x)\phi^{(2)}(x)\,dx \nonumber\\ &+ \int_{x_1}^{x_2} b(t, x) u_x(t, x)\phi^{(2)}(x)\,dx \nonumber\\ &+ \int_{x_1}^{x_2} c(t, x) u(t, x)\phi^{(2)}(x)\,dx,\nonumber\\ \int_{x_2}^{x_3} u_t(t, x)\phi^{(3)}(x)\,dx =& \int_{x_2}^{x_3} a(t, x)u_{xx}(t, x)\phi^{(3)}(x)\,dx \nonumber\\ &+ \int_{x_2}^{x_3} b(t, x) u_x(t, x)\phi^{(3)}(x)\,dx \nonumber\\ &+ \int_{x_2}^{x_3} c(t, x) u(t, x)\phi^{(3)}(x)\,dx,\nonumber \end{align*} where we write \begin{align*} \phi(x) = \begin{cases} \phi^{(1)}(x) \quad x\in[x_0, x_1]\\ \phi^{(2)}(x) \quad x\in[x_1, x_2]\\ \phi^{(3)}(x) \quad x\in[x_2, x_3]\\ \end{cases}. \end{align*} Note that at the interface we have \(\phi^{(1)}(x_1) = \phi^{(2)}(x_1)\), \(\phi^{(2)}(x_2) = \phi^{(3)}(x_2)\).

\begin{align*} u(t, \xi) \approx& \sum_{p=0}^{N_1} U_p^{(1)}(t)\ell_p^{(1)}(\xi^{(1)}),\\ \phi^{(1)}(\xi) \approx& \sum_{i=0}^{N_1} \phi_i^{(1)}\ell_i^{(1)}(\xi^{(1)}),\\ \int_{-1}^1 f(\xi)\,d\xi \approx& \sum_{n=0}^{N_1}f(\xi_n)w_n^{(1)}, \end{align*}

\begin{align*} \frac{d_1}2~ \sum_{i=0}^{N_1} \dot U_i^{(1)}(t)\phi_i^{(1)}w_i^{(1)} =& -\frac2{d_1}~ \sum_{p=0}^{N_1}\sum_{i=0}^{N_1}\sum_{n=0}^{N_1} U_p^{(1)}D_{np}^{(1)}A_i^{(1)}\phi_i^{(1)}D_{ni}^{(1)}w_n^{(1)} \nonumber\\ &+~ \sum_{p=0}^{N_1}\sum_{i=0}^{N_1} U_p^{(1)}D_{ip}^{(1)}(b^{(1)}-a_x^{(1)})_i\phi_i^{(1)}w_i^{(1)} \nonumber\\ &+~ \frac{d_1}2~\sum_{i=0}^{N_1} U_i^{(1)}C_i^{(1)}\phi_i^{(1)}w_i^{(1)}, \label{eq:phiieqelem1} \end{align*}

\begin{align*} &\frac{d_1}2~ \sum_{i=0}^{N_1} \dot U_i^{(1)}(t)\phi_i^{(1)}w_i^{(1)} + \frac{d_2}2~ \sum_{i=0}^{N_2} \dot U_i^{(2)}(t)\phi_i^{(2)}w_i^{(2)} + \frac{d_3}2~ \sum_{i=0}^{N_3} \dot U_i^{(3)}(t)\phi_i^{(3)}w_i^{(3)}\nonumber\\ &= \left(-\frac2{d_1}~ \sum_{p=0}^{N_1}\sum_{i=0}^{N_1}\sum_{n=0}^{N_1} U_p^{(1)}D_{np}^{(1)}A_i^{(1)}\phi_i^{(1)}D_{ni}^{(1)}w_n^{(1)} \right.\nonumber\\ &\qquad\left. +~ \sum_{p=0}^{N_1}\sum_{i=0}^{N_1} U_p^{(1)}D_{ip}^{(1)}(b^{(1)}-a_x^{(1)})_i\phi_i^{(1)}w_i^{(1)} +~ \frac{d_1}2~\sum_{i=0}^{N_1} U_i^{(1)}C_i^{(1)}\phi_i^{(1)}w_i^{(1)}\right)\nonumber\\ &\quad+ \left(-\frac2{d_2}~ \sum_{p=0}^{N_2}\sum_{i=0}^{N_2}\sum_{n=0}^{N_2} U_p^{(2)}D_{np}^{(2)}A_i^{(2)}\phi_i^{(2)}D_{ni}^{(2)}w_n^{(2)} \right.\nonumber\\ &\qquad\quad\left. +~ \sum_{p=0}^{N_2}\sum_{i=0}^{N_2} U_p^{(2)}D_{ip}^{(2)}(b^{(2)}-a_x^{(2)})_i\phi_i^{(2)}w_i^{(2)} +~ \frac{d_2}2~\sum_{i=0}^{N_2} U_i^{(2)}C_i^{(2)}\phi_i^{(2)}w_i^{(2)}\right)\nonumber\\ &\quad+ \left(-\frac2{d_3}~ \sum_{p=0}^{N_3}\sum_{i=0}^{N_3}\sum_{n=0}^{N_3} U_p^{(3)}D_{np}^{(3)}A_i^{(3)}\phi_i^{(3)}D_{ni}^{(3)}w_n^{(3)} \right.\nonumber\\ &\qquad\quad\left. +~ \sum_{p=0}^{N_3}\sum_{i=0}^{N_3} U_p^{(3)}D_{ip}^{(3)}(b^{(3)}-a_x^{(3)})_i\phi_i^{(3)}w_i^{(3)} +~\frac{d_3}2~\sum_{i=0}^{N_3} U_i^{(3)}C_i^{(3)}\phi_i^{(3)}w_i^{(3)}\right).\qquad\qquad\label{eq:fullsem} \end{align*} and the identities at the interface translate to \(\phi_{N_1}^{(1)} = \phi_{0}^{(2)}\), \(\phi_{N_2}^{(2)} = \phi_{0}^{(3)}\).

Setting \(\phi_{i^*}^{(1)} = 1\) for a specific \(i^*\) and all other \(\phi_{i}^{(j)}\)’s to zero: \begin{align*} \frac{d_1}{2}~\dot U_{i^*}^{(1)}(t) w_{i^*}^{(1)} =& -\frac2{d_1}~ \sum_{p=0}^{N_1}\sum_{n=0}^{N_1} U_p^{(1)}(t)D_{np}^{(1)}A_{i^*}^{(1)}D_{n{i^*}}^{(1)}w_n^{(1)} \nonumber\\ &+~ \sum_{p=0}^{N_1} U_p^{(1)}(t)D_{{i^*}p}^{(1)}(b-a_x)_{i^*}^{(1)}w_{i^*}^{(1)}\nonumber\\ &+~ \frac{d_1}2~U_{i^*}^{(1)}(t)C_{i^*}^{(1)}w_{i^*}^{(1)}.\label{eq:ODEe1} \end{align*} Similarly \begin{align*} \frac{d_2}{2}~\dot U_{i}^{(2)}(t) w_{i}^{(2)} =& -\frac2{d_2}~ \sum_{p=0}^{N_2}\sum_{n=0}^{N_2} U_p^{(2)}(t)D_{np}^{(2)}A_{i}^{(2)}D_{n{i}}^{(2)}w_n^{(2)} \nonumber\\ &+~ \sum_{p=0}^{N_2} U_p^{(2)}(t)D_{{i}p}^{(2)}(b-a_x)_{i}^{(2)}w_{i}^{(2)}\nonumber\\ &+~ \frac{d_2}2~ U_{i}^{(2)}(t)C_{i}^{(2)}w_{i}^{(2)},\qquad\forall i=1, 2, \ldots, N_2-1, \label{eq:ODEe2}\\ \frac{d_3}{2}~\dot U_{i}^{(3)}(t) w_{i}^{(3)} =& -\frac2{d_3}~ \sum_{p=0}^{N_3}\sum_{n=0}^{N_3} U_p^{(3)}(t)D_{np}^{(3)}A_{i}^{(3)}D_{n{i}}^{(3)}w_n^{(3)} \nonumber\\ &+~ \sum_{p=0}^{N_3} U_p^{(3)}(t)D_{{i}p}^{(3)}(b-a_x)_{i}^{(3)}w_{i}^{(3)}\nonumber\\ &+~ \frac{d_3}2~U_{i}^{(3)}(t)C_{i}^{(3)}w_{i}^{(3)}, \qquad\forall i=1, 2, \ldots, N_3-1.\label{eq:ODEe3} \end{align*}

Setting \(\phi_{N_1}^{(1)} = \phi_{0}^{(2)} = 1\), all other \(\phi_{i}^{(j)}\)’s to zero: \begin{align*} &\frac{d_1}2~ \dot U_{N_1}^{(1)}(t)w_{N_1}^{(1)} + \frac{d_2}2~ \dot U_0^{(2)}(t)w_0^{(2)}\nonumber\\ &= \left(-\frac2{d_1}~ \sum_{p=0}^{N_1}\sum_{n=0}^{N_1} U_p^{(1)}D_{np}^{(1)}A_{N_1}^{(1)}D_{n{N_1}}^{(1)}w_n^{(1)} \right.\nonumber\\ &\qquad\left. +~ \sum_{p=0}^{N_1} U_p^{(1)}D_{{N_1}p}^{(1)}(b^{(1)}-a_x^{(1)})_{N_1}w_{N_1}^{(1)} +~ \frac{d_1}2~ U_{N_1}^{(1)}C_{N_1}^{(1)}w_{N_1}^{(1)}\right)\nonumber\\ &\quad+ \left(-\frac2{d_2}~ \sum_{p=0}^{N_2}\sum_{n=0}^{N_2} U_p^{(2)}D_{np}^{(2)}A_0^{(2)}D_{n0}^{(2)}w_n^{(2)} \right.\nonumber\\ &\qquad\quad\left. +~ \sum_{p=0}^{N_2} U_p^{(2)}D_{0p}^{(2)}(b^{(2)}-a_x^{(2)})_0\phi_0^{(2)}w_0^{(2)} +~ \frac{d_2}2~ U_0^{(2)}C_0^{(2)}\phi_0^{(2)}w_0^{(2)}\right),\label{eq:ODEx1} \end{align*}

Setting \(\phi_{N_2}^{(2)} = \phi_{0}^{(3)} = 1\), all other \(\phi_{i}^{(j)}\)’s to zero: \begin{align*} &\frac{d_2}2~ \dot U_{N_2}^{(2)}(t)w_{N_2}^{(2)} + \frac{d_3}2~ \dot U_0^{(3)}(t)w_0^{(3)}\nonumber\\ &= \left(-\frac2{d_2}~ \sum_{p=0}^{N_2}\sum_{n=0}^{N_2} U_p^{(2)}D_{np}^{(2)}A_{N_2}^{(2)}D_{n{N_2}}^{(2)}w_n^{(2)} \right.\nonumber\\ &\qquad\left. +~ \sum_{p=0}^{N_2} U_p^{(2)}D_{{N_2}p}^{(2)}(b^{(2)}-a_x^{(2)})_{N_2}w_{N_2}^{(2)} +~ \frac{d_2}2~ U_{N_2}^{(2)}C_{N_2}^{(2)}w_{N_2}^{(2)}\right)\nonumber\\ &\quad+ \left(-\frac2{d_3}~ \sum_{p=0}^{N_3}\sum_{n=0}^{N_3} U_p^{(3)}D_{np}^{(3)}A_0^{(3)}D_{n0}^{(3)}w_n^{(3)} \right.\nonumber\\ &\qquad\quad\left. +~ \sum_{p=0}^{N_3} U_p^{(3)}D_{0p}^{(3)}(b^{(3)}-a_x^{(3)})_0\phi_0^{(3)}w_0^{(3)} +~ \frac{d_3}2~ U_0^{(3)}C_0^{(3)}\phi_0^{(3)}w_0^{(3)}\right),\label{eq:ODEx2} \end{align*}

\begin{align*} \Omega \dot{\vec U}(t) = M \vec U(t), \end{align*}

(Pictures that doesn’t display on Jupyter) where the overlapping region is the sum of two elements form two matrices, For example the overlapping region of \(M^{(1)}\) and \(M^{(2)}\) is \(M^{(1)}_{N_1 N_1} + M^{(2)}_{0 0}\). And

\begin{align*} M^{(i)} = \left( -\frac2{d_i} (D^{(i)})^TA^{(i)}W^{(i)}D^{(i)} + \left(B^{(i)}-A_x^{(i)}\right)W^{(i)}D^{(i)} + \frac{d_i}2C^{(i)}W^{(i)}D^{(i)}\right), \end{align*}

:nbsphinx-math:`begin{align*} &A_x^{(i)} = left( begin{array}{ccccc} a_x(t, xi_0^{(i)}) & & & & 0 \ & a_x(t, xi_1^{(i)})& & & \ & & a_x(t, xi_2^{(i)})& & \ & & & ddots & \ 0 & & & & a_x(t, xi_{N_i}^{(i)}) end{array} right), \\ &A^{(i)} = left( begin{array}{ccccc} a(t, xi_0^{(i)}) & & & & 0 \ & a(t, xi_1^{(i)})& & & \ & & a(t, xi_2^{(i)})& & \ & & & ddots & \ 0 & & & & a(t, xi_{N_i}^{(i)}) end{array} right), B^{(i)} = left( begin{array}{ccccc} b(t, xi_0^{(i)}) & & & & 0 \ & b(t, xi_1^{(i)})& & & \ & & b(t, xi_2^{(i)})& & \ & & & ddots & \ 0 & & & & b(t, xi_{N_i}^{(i)}) end{array} right), \\ &C^{(i)} = left( begin{array}{ccccc} c(t, xi_0^{(i)}) & & & & 0 \ & c(t, xi_1^{(i)})& & & \ & & c(t, xi_2^{(i)})& & \ & & & ddots & \ 0 & & & & c(t, xi_{N_i}^{(i)}) end{array} right), W^{(i)} = left( begin{array}{ccccc} w_0^{(i)} & & & & 0 \ & w_1^{(i)} & & & \ & & w_2^{(i)} & & \ & & & ddots & \ 0 & & & & w_{N_i}^{(i)} end{array} right), \\ &D^{(i)} = left( begin{array}{ccccc}

D_{00}^{(i)} & D_{01}^{(i)} & D_{02}^{(i)} & cdots & D_{0{N_i}}^{(i)} \ D_{10}^{(i)} & D_{11}^{(i)} & D_{12}^{(i)} & cdots & D_{1{N_i}}^{(i)} \ D_{20}^{(i)} & D_{21}^{(i)} & D_{22}^{(i)} & cdots & D_{2{N_i}}^{(i)} \

vdots & vdots & vdots & ddots & vdots \

D_{{N_i}0}^{(i)} & D_{{N_i}1}^{(i)} & D_{{N_i}2}^{(i)} & cdots & D_{{N_i}{N_i}}^{(i)}

end{array} right) = left( begin{array}{ccccc}

frac{d}{dxi^{(i)}}ell_0^{(i)}(xi_0^{(i)}) & frac{d}{dxi^{(i)}}ell_1^{(i)}(xi_0^{(i)}) & frac{d}{dxi^{(i)}}ell_2^{(i)}(xi_0^{(i)}) & cdots & frac{d}{dxi^{(i)}}ell_{N_i}^{(i)}(xi_0^{(i)}) \ frac{d}{dxi^{(i)}}ell_0^{(i)}(xi_1^{(i)}) & frac{d}{dxi^{(i)}}ell_1^{(i)}(xi_1^{(i)}) & frac{d}{dxi^{(i)}}ell_2^{(i)}(xi_1^{(i)}) & cdots & frac{d}{dxi^{(i)}}ell_{N_i}^{(i)}(xi_1^{(i)}) \ frac{d}{dxi^{(i)}}ell_0^{(i)}(xi_2^{(i)}) & frac{d}{dxi^{(i)}}ell_1^{(i)}(xi_2^{(i)}) & frac{d}{dxi^{(i)}}ell_2^{(i)}(xi_2^{(i)}) & cdots & frac{d}{dxi^{(i)}}ell_{N_i}^{(i)}(xi_2^{(i)}) \

vdots & vdots & vdots & ddots & vdots \

frac{d}{dxi^{(i)}}ell_0^{(i)}(xi_{N_i}^{(i)}) & frac{d}{dxi^{(i)}}ell_1^{(i)}(xi_{N_i}^{(i)}) & frac{d}{dxi^{(i)}}ell_2^{(i)}(xi_{N_i}^{(i)}) & cdots & frac{d}{dxi^{(i)}}ell_{N_i}^{(i)}(xi_{N_i}^{(i)})

end{array} right). end{align*}`

\begin{align*} \vec U(t_{n+1}) &= \vec U(t_n) + \int_{t_n}^{t_{n+1}} \dot{\vec U}(t)\,dt\\ &= \vec U(t_n) + \int_{t_n}^{t_{n+1}} L(t)\vec U(t)\,dt. \end{align*}

\begin{align*} \vec U(t_{n+1}) = \vec U(t_n) + \frac{\Delta t}{2}\left( L(t_n)\vec U(t_n) + L(t_{n+1})\vec U(t_{n+1}) \right). \end{align*} Rearrange the terms to obtain \begin{align*} \left(I - \frac{\Delta t}{2}L(t_{n+1})\right)\vec U(t_{n+1}) = \left(I + \frac{\Delta t}{2}L(t_{n})\right)\vec U(t_n). \end{align*}