CFD Examples 1D

Runge-Kutta

Consider an ordinary differential equation of the form \(\cfrac{dy}{dx} = f(x, y)\) with initial condition \(y(x_0) = y_0\). For this, we can define the formulas for Runge-Kutta methods as follows.

Third-order version1:

\[\begin{split}\begin{array}{l} k_1=f(x_{n},y_{n})\\ k_2=f(x_{n}+\cfrac{1}{2}h,y_{n}+\cfrac{1}{2}hk_{1})\\ k_3=f(x_{0}+h,y_{n}-hk_{1}+2hk_{2})\\ y_{n+1}=y_{n}+\cfrac{1}{6}h(k_{1}+4k_{2}+k_{3})\\ \end{array}\end{split}\]

Third-order version2:

\[\begin{split}\begin{array}{l} k_1=f(x_{n},y_{n})\\ k_2=f(x_{n}+\cfrac{1}{3}h,y_{n}+\cfrac{1}{3}hk_{1})\\ k_3=f(x_{n}+\cfrac{2}{3}h,y_{n}+\cfrac{2}{3}hk_{2})\\ y_{n+1}=y_{n}+\cfrac{1}{4}h(k_{1}+3k_{3})\\ \end{array}\end{split}\]

Third-order version3: (Ralston formula)

\[\begin{split}\begin{array}{l} k_1=f(x_{n},y_{n})\\ k_2=f(x_{n}+\cfrac{1}{2}h,y_{n}+\cfrac{1}{2}hk_{1})\\ k_3=f(x_{n}+\cfrac{3}{4}h,y_{n}+\cfrac{3}{4}hk_{2})\\ y_{n+1}=y_{n}+\cfrac{1}{9}h(2k_{1}+3k_{2}+4k_{3})\\ \end{array}\end{split}\]

Third order TVD Runge-Kutta method

\[\begin{split}\begin{array}{l} u^{(1)}=u^{n}+\Delta tL(u^{n})\\ u^{(2)}=\cfrac{3}{4}u^{n}+\cfrac{1}{4}u^{(1)}+\cfrac{1}{4}\Delta tL(u^{(1)})\\ u^{n+1}=\cfrac{1}{3}u^{n}+\cfrac{2}{3}u^{(2)}+\cfrac{2}{3}\Delta tL(u^{(2)})\\ \end{array}\end{split}\]

Let us split the interval \([x_{0}, x_{*}]\) into \(n\) equal segments of length \(h = (x_{*} − x_{0})/n\). Let \(y_{1}, . . . , y_{n}\) denote the approximate values of the function \(y(x)\) at the partitioning points \(x_{1}, . . . , x_{n} = x_{*}\). The discrete set of points \(X_{h} = \{ x_{0}, x_{1}, \dots, x_{n}\}\) is called a mesh, the individual points \(x_{k}\) are mesh nodes, and \(h\) is a mesh increment or step size.

\[x_{k} = x_{0} + kh\]

Heat Equation

\[\begin{split}\begin{array}{l} \cfrac{\partial u}{ \partial t}=\alpha \cfrac{\partial ^{2}u}{ \partial x^{2}}=f(t,u)=g(u)\\ u\sim y\rightarrow \cfrac{\partial y}{ \partial t}=\alpha \cfrac{\partial ^{2}y}{ \partial x^{2}}=f(t,y)=g(y)\\ \end{array}\end{split}\]

Third-order Runge-Kutta version1:

\[\begin{split}\begin{array}{l} h=\Delta t\\ \displaystyle k_{1}=f(t_{n},u_{n})=\alpha \cfrac{\partial ^{2}u}{\partial x^{2}}\Bigg|_{(t_{n},u_{n})} =L(u_{n})\\ \displaystyle k_{2}=f(t_{n}+\frac{1}{2}h,u_{n}+\frac{1}{2}hk_{1})=\alpha \cfrac{\partial ^{2}u}{\partial x^{2}}\Bigg|_{(t_{n}+\frac{1}{2}h,u_{n}+\frac{1}{2}hk_{1})} =L(u_{n}+\frac{1}{2}hk_{1})\\ \displaystyle k_{3}=f(t_{n}+h,u_{n}+h(-k_{1}+2k_{2}))=\alpha \cfrac{\partial ^{2}u}{\partial x^{2}}\Bigg|_{(t_{n}+h,u_{n}+h(-k_{1}+2k_{2}))} =L(u_{n}+h(-k_{1}+2k_{2}))\\ \displaystyle u_{n+1}=u_{n}+\frac{1}{6}h(k_{1}+4k_{2}+k_{3}) \end{array}\end{split}\]

Specifically

\[\begin{split}\begin{array}{l} h=\Delta t\\ \displaystyle k_{1}=\alpha \cfrac{u^{n}_{i+1}-2u^{n}_{i}+u^{n}_{i-1}}{\Delta x^{2}}=L(u_{n})\\ u^{(1)}_i=u_{n}+\frac{1}{2}hk_{1}=u_{n}+\frac{1}{2}\Delta tL(u_{n})\\ \displaystyle k_{2}=\alpha \cfrac{u^{(1)}_{i+1}-2u^{(1)}_{i}+u^{(1)}_{i-1}}{\Delta x^{2}}=L(u_{n}+\frac{1}{2}hk_{1})=L(u^{(1)})\\ u^{(2)}_i=u_{n}+h(-k_{1}+2k_{2})\\ \displaystyle k_{3}=\alpha \cfrac{u^{(2)}_{i+1}-2u^{(2)}_{i}+u^{(2)}_{i-1}}{\Delta x^{2}} =L(u^{(2)})\\ \displaystyle u_{n+1}=u_{n}+\frac{1}{6}h(k_{1}+4k_{2}+k_{3}) \end{array}\end{split}\]

\[\begin{split}\begin{array}{l} h=\Delta t\\ \displaystyle \hat{k}_{1}=hk_{1}=\alpha \Delta t\cfrac{u^{n}_{i+1}-2u^{n}_{i}+u^{n}_{i-1}}{\Delta x^{2}}=\Delta tL(u_{n})\\ u^{(1)}_i=u_{n}+\frac{1}{2}\hat{k}_{1}=u_{n}+\frac{1}{2}\Delta tL(u_{n})\\ \displaystyle \hat{k}_{2}=hk_{2}=\alpha \Delta t\cfrac{u^{(1)}_{i+1}-2u^{(1)}_{i}+u^{(1)}_{i-1}}{\Delta x^{2}}=\Delta tL(u_{n}+\frac{1}{2}hk_{1})=\Delta tL(u^{(1)})\\ u^{(2)}_i=u_{n}+h(-k_{1}+2k_{2})=u_{n}+(-\hat{k}_{1}+2\hat{k}_{2})\\ \displaystyle \hat{k}_{3}=hk_{3}=\alpha \Delta t\cfrac{u^{(2)}_{i+1}-2u^{(2)}_{i}+u^{(2)}_{i-1}}{\Delta x^{2}} =\Delta tL(u^{(2)})\\ \displaystyle u_{n+1}=u_{n}+\frac{1}{6}(\hat{k}_{1}+4\hat{k}_{2}+\hat{k}_{3}) \end{array}\end{split}\]

Third-order Runge-Kutta version2:

\[\begin{split}\begin{array}{l} h=\Delta t\\ \displaystyle \hat{k}_{1}=hf(t_{n},u_{n})=\alpha \Delta t\cfrac{\partial ^{2}u}{\partial x^{2}}\Bigg|_{(t_{n},u_{n})} =\Delta tL(u_{n})\\ \displaystyle \hat{k}_{2}=hf(t_{n}+\frac{1}{3}h,u_{n}+\frac{1}{3}\hat{k}_{1})=\alpha\Delta t\cfrac{\partial ^{2}u}{\partial x^{2}}\Bigg|_{(t_{n}+\frac{1}{3}h,u_{n}+\frac{1}{3}\hat{k}_{1})} =\Delta tL(u_{n}+\frac{1}{3}\hat{k}_{1})\\ \displaystyle \hat{k}_{3}=hf(t_{n}+\frac{2}{3}h,u_{n}+\frac{2}{3}\hat{k}_{2})=\alpha\Delta t\cfrac{\partial ^{2}u}{\partial x^{2}}\Bigg|_{(t_{n}+\frac{2}{3}h,u_{n}+\frac{2}{3}\hat{k}_{2})} =\Delta tL(u_{n}+\frac{2}{3}\hat{k}_{2})\\ \displaystyle u_{n+1}=u_{n}+\frac{1}{4}(\hat{k}_{1}+3\hat{k}_{3}) \end{array}\end{split}\]

\[\begin{split}\begin{array}{l} h=\Delta t\\ \displaystyle \hat{k}_{1}=hk_{1}=\alpha \Delta t\cfrac{u^{n}_{i+1}-2u^{n}_{i}+u^{n}_{i-1}}{\Delta x^{2}}=\Delta tL(u_{n})\\ u^{(1)}_i=u_{n}+\frac{1}{3}\hat{k}_{1}=u_{n}+\frac{1}{3}\Delta tL(u_{n})\\ \displaystyle \hat{k}_{2}=hk_{2}=\alpha \Delta t\cfrac{u^{(1)}_{i+1}-2u^{(1)}_{i}+u^{(1)}_{i-1}}{\Delta x^{2}}=\Delta tL(u_{n}+\frac{1}{3}hk_{1})=\Delta tL(u^{(1)})\\ u^{(2)}_i=u_{n}+\frac{2}{3}\hat{k}_{2}\\ \displaystyle \hat{k}_{3}=hk_{3}=\alpha \Delta t\cfrac{u^{(2)}_{i+1}-2u^{(2)}_{i}+u^{(2)}_{i-1}}{\Delta x^{2}} =\Delta tL(u^{(2)})\\ \displaystyle u_{n+1}=u_{n}+\frac{1}{4}(\hat{k}_{1}+3\hat{k}_{3}) \end{array}\end{split}\]

Third-order Runge-Kutta version3(Ralston formula):

\[\begin{split}\begin{array}{l} h=\Delta t\\ \displaystyle \hat{k}_{1}=hf(t_{n},u_{n})=\alpha \Delta t\cfrac{\partial ^{2}u}{\partial x^{2}}\Bigg|_{(t_{n},u_{n})} =\Delta tL(u_{n})\\ \displaystyle \hat{k}_{2}=hf(t_{n}+\frac{1}{2}h,u_{n}+\frac{1}{2}\hat{k}_{1})=\alpha\Delta t\cfrac{\partial ^{2}u}{\partial x^{2}}\Bigg|_{(t_{n}+\frac{1}{2}h,u_{n}+\frac{1}{2}\hat{k}_{1})} =\Delta tL(u_{n}+\frac{1}{2}\hat{k}_{1})\\ \displaystyle \hat{k}_{3}=hf(t_{n}+\frac{3}{4}h,u_{n}+\frac{3}{4}\hat{k}_{2})=\alpha\Delta t\cfrac{\partial ^{2}u}{\partial x^{2}}\Bigg|_{(t_{n}+\frac{3}{4}h,u_{n}+\frac{3}{4}\hat{k}_{2})} =\Delta tL(u_{n}+\frac{3}{4}\hat{k}_{2})\\ \displaystyle u_{n+1}=u_{n}+\frac{1}{9}(2\hat{k}_{1}+3\hat{k}_{2}+4\hat{k}_{3}) \end{array}\end{split}\]

\[\begin{split}\begin{array}{l} h=\Delta t\\ \displaystyle \hat{k}_{1}=hk_{1}=\alpha \Delta t\cfrac{u^{n}_{i+1}-2u^{n}_{i}+u^{n}_{i-1}}{\Delta x^{2}}=\Delta tL(u_{n})\\ u^{(1)}_i=u_{n}+\frac{1}{2}\hat{k}_{1}=u_{n}+\frac{1}{2}\Delta tL(u_{n})\\ \displaystyle \hat{k}_{2}=hk_{2}=\alpha \Delta t\cfrac{u^{(1)}_{i+1}-2u^{(1)}_{i}+u^{(1)}_{i-1}}{\Delta x^{2}}=\Delta tL(u_{n}+\frac{1}{2}hk_{1})=\Delta tL(u^{(1)})\\ u^{(2)}_i=u_{n}+\frac{3}{4}\hat{k}_{2}\\ \displaystyle \hat{k}_{3}=hk_{3}=\alpha \Delta t\cfrac{u^{(2)}_{i+1}-2u^{(2)}_{i}+u^{(2)}_{i-1}}{\Delta x^{2}} =\Delta tL(u^{(2)})\\ \displaystyle u_{n+1}=u_{n}+\frac{1}{9}(2\hat{k}_{1}+3\hat{k}_{2}+4\hat{k}_{3}) \end{array}\end{split}\]

Third order TVD Runge-Kutta method

\[\begin{split}\begin{array}{l} \hat{k}_{1}=\alpha\Delta t\cfrac{u_{i+1}^{(n)}-2 u_{i}^{(n)}+u_{i-1}^{(n)}}{\Delta x^{2}}\\ u_{i}^{(1)} = u_{i}^{(n)}+\hat{k}_{1}\\ \hat{k}_{2}=\alpha\Delta t\cfrac{u_{i+1}^{(1)}-2 u_{i}^{(1)}+u_{i-1}^{(1)}}{\Delta x^{2}}\\ u_{i}^{(2)} = u_{i}^{(n)}+\frac{1}{4} \hat{k}_{1}+\frac{1}{4}\hat{k}_{2} \\ \hat{k}_{3}=\alpha\Delta t\cfrac{u_{i+1}^{(2)}-2 u_{i}^{(2)}+u_{i-1}^{(2)}}{\Delta x^{2}}\\ u_{i}^{(n+1)} = u_{i}^{(n)}+\frac{1}{6} (\hat{k}_{1}+\hat{k}_{2}+4\hat{k}_{3})\ \end{array}\end{split}\]

Equivalent to

\[\begin{split}\begin{align} u_{i}^{(1)} & = u_{i}^{(n)}+\frac{\alpha \Delta t}{\Delta x^{2}}\left(u_{i+1}^{(n)}-2 u_{i}^{(n)}+u_{i-1}^{(n)}\right),\\ u_{i}^{(2)} & = \frac{3}{4} u_{i}^{(n)}+\frac{1}{4} u_{i}^{(1)}+\frac{1}{4}\frac{\alpha \Delta t}{ \Delta x^{2}}\left(u_{i+1}^{(1)}-2 u_{i}^{(1)}+u_{i-1}^{(1)}\right) \\ u_{i}^{(n+1)} & = \frac{1}{3} u_{i}^{(n)}+\frac{2}{3} u_{i}^{(2)}+\frac{2}{3}\frac{\alpha \Delta t}{\Delta x^{2}}\left(u_{i+1}^{(2)}-2 u_{i}^{(2)}+u_{i-1}^{(2)}\right) . \end{align}\end{split}\]

\[\begin{split}\begin{array}{l} u^{(1)}=u^{n}+\Delta tL(u^{n})\\ u^{(2)}=\cfrac{3}{4}u^{n}+\cfrac{1}{4}u^{(1)}+\cfrac{1}{4}\Delta tL(u^{(1)})\\ u^{n+1}=\cfrac{1}{3}u^{n}+\cfrac{2}{3}u^{(2)}+\cfrac{2}{3}\Delta tL(u^{(2)})\\ \end{array}\end{split}\]

Heat Equation

The one-dimensional heat equation is given as

\[\cfrac{\partial u}{ \partial t}=\alpha \cfrac{\partial ^{2}u}{ \partial x^{2}}\]

where \(u\) is the field variable, \(t\) is the time variable, and \(\alpha\) is the diffusivity of the medium. The heat equation describes the evolution of the field variable over time in a medium.

Forward Time Central Space (FTCS) Scheme

\[\frac{u_{i}^{(n+1)}-u_{i}^{(n)}}{\Delta t}=\alpha \frac{u_{i+1}^{(n)}-2 u_{i}^{(n)}+u_{i-1}^{(n)}}{\Delta x^{2}},\]

we can re-write the above equation as an explicit update formula

\[u_{i}^{(n+1)}=u_{i}^{(n)}+\alpha \cfrac{\Delta t}{\Delta x^{2}}({u_{i+1}^{(n)}-2 u_{i}^{(n)}+u_{i-1}^{(n)}})\]

We use the computational domain \(x\in [-1,1]\) and \(\alpha=1/\pi^{2}\). The initial condition is \(u(t=0,x)=-sin(\pi x)\). The analytical solution to the one-dimensional heat equation is given by

\[u(t,x)=-e^{-t}sin(\pi x)\]

We use \(\Delta x=0.025\) and \(\Delta 𝑡=0.0025\) for spatial and temporal discretization.

Runge-Kutta Numerical Scheme

Runge-Kutta methods tries to improve the accuracy of temporal term by evaluating \(f\) at intermediate points between \(t_n\) and \(t_{n+1}\). The additional steps lead to an increase in computational time, but the temporal accuracy is increased. The time integration of the heat equation using third-order Runge-Kutta scheme is given below:

\[\begin{split}\begin{align} u_{i}^{(1)} & = u_{i}^{(n)}+\frac{\alpha \Delta t}{\Delta x^{2}}\left(u_{i+1}^{(n)}-2 u_{i}^{(n)}+u_{i-1}^{(n)}\right),\\ u_{i}^{(2)} & = \frac{3}{4} u_{i}^{(n)}+\frac{1}{4} u_{i}^{(1)}+\frac{1}{4}\frac{\alpha \Delta t}{ \Delta x^{2}}\left(u_{i+1}^{(1)}-2 u_{i}^{(1)}+u_{i-1}^{(1)}\right) \\ u_{i}^{(n+1)} & = \frac{1}{3} u_{i}^{(n)}+\frac{2}{3} u_{i}^{(2)}+\frac{2}{3}\frac{\alpha \Delta t}{\Delta x^{2}}\left(u_{i+1}^{(2)}-2 u_{i}^{(2)}+u_{i-1}^{(2)}\right) . \end{align}\end{split}\]

Crank–Nicolson method

Crank–Nicolson method

The Crank–Nicolson method is based on the trapezoidal rule, giving second-order convergence in time. For example, in one dimension, suppose the partial differential equation is

\[\frac{\partial u}{\partial t}=F\left(u, x, t, \frac{\partial u}{\partial x}, \frac{\partial^{2} u}{\partial x^{2}}\right) .\]

Letting \(u(i\Delta x, n\Delta t)=u^{n}_{i}\) and \(F^{n}_{i}=F\) evaluated for \(i\), \(n\) and \(u^{n}_{i}\), the equation for Crank-Nicolson method is a combination of the forward Euler method at \(n\) and the backward Euler method at \(n+1\)

\[\begin{split}\begin{array}{|l|l|} \hline \frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t}=F_{i}^{n}\left(u, x, t, \frac{\partial u}{\partial x}, \frac{\partial^{2} u}{\partial x^{2}}\right) & \text { forward Euler } \\ \hline \frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t}=F_{i}^{n+1}\left(u, x, t, \frac{\partial u}{\partial x}, \frac{\partial^{2} u}{\partial x^{2}}\right) & \text { backward Euler } \\ \hline \frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t}=\frac{1}{2}\left[F_{i}^{n+1}\left(u, x, t, \frac{\partial u}{\partial x}, \frac{\partial^{2} u}{\partial x^{2}}\right)+F_{i}^{n}\left(u, x, t, \frac{\partial u}{\partial x}, \frac{\partial^{2} u}{\partial x^{2}}\right)\right] & \text { Crank-Nicolson } \\ \hline \end{array}\end{split}\]

Note that this is an implicit method: to get the “next” value of u in time, a system of algebraic equations must be solved. If the partial differential equation is nonlinear, the discretization will also be nonlinear, so that advancing in time will involve the solution of a system of nonlinear algebraic equations, though linearizations are possible. In many problems, especially linear diffusion, the algebraic problem is tridiagonal and may be efficiently solved with the tridiagonal matrix algorithm.

Example: 1D diffusion

The Crank–Nicolson method is often applied to diffusion problems. As an example, for linear diffusion,

\[\cfrac{\partial u}{ \partial t}=\alpha \cfrac{\partial ^{2}u}{ \partial x^{2}}\]

applying a finite difference spatial discretization for the right-hand side, the Crank–Nicolson discretization is then

\[\cfrac{u^{n+1}_{i}-u^{n}_{i}}{\Delta t}=\cfrac{1}{2}\cfrac{\alpha}{\Delta x^{2}} \left[(u^{n+1}_{i+1}-2u^{n+1}_{i}+u^{n+1}_{i-1})+(u^{n}_{i+1}-2u^{n}_{i}+u^{n}_{i-1})\right]\]

or, letting \(r=\cfrac{1}{2}\cfrac{\alpha\Delta t}{\Delta x^{2}}\),

\[-ru^{n+1}_{i+1}+(1+2r)u^{n+1}_{i}-ru^{n+1}_{i-1}=ru^{n}_{i+1}+(1-2r)u^{n}_{i}+ru^{n}_{i-1}\]

\[-ru^{n+1}_{i-1}+(1+2r)u^{n+1}_{i}-ru^{n+1}_{i+1}=ru^{n}_{i-1}+(1-2r)u^{n}_{i}+ru^{n}_{i+1}\]

\[\begin{split}\begin{array}{l} 1:-ru^{n+1}_{0}+(1+2r)u^{n+1}_{1}-ru^{n+1}_{2}=d_{1}\\ 2:-ru^{n+1}_{1}+(1+2r)u^{n+1}_{2}-ru^{n+1}_{3}=d_{2}\\ 3:-ru^{n+1}_{2}+(1+2r)u^{n+1}_{3}-ru^{n+1}_{4}=d_{3}\\ \cdots \\ i:-ru^{n+1}_{i-1}+(1+2r)u^{n+1}_{i}-ru^{n+1}_{i+1}=d_{i}\\ \cdots \\ N-2:-ru^{n+1}_{N-3}+(1+2r)u^{n+1}_{N-2}-ru^{n+1}_{N-1}=d_{N-2}\\ N-1:-ru^{n+1}_{N-2}+(1+2r)u^{n+1}_{N-1}-ru^{n+1}_{N}=d_{N-1}\\ \end{array}\end{split}\]

\[d_{i}=ru^{n}_{i-1}+(1-2r)u^{n}_{i}+ru^{n}_{i+1}\]

We can rewrite this as

\[\begin{split}\begin{array}{l} 1:(1+2r)u^{n+1}_{1}-ru^{n+1}_{2}=d_{1}+ru^{n+1}_{0}\\ 2:-ru^{n+1}_{1}+(1+2r)u^{n+1}_{2}-ru^{n+1}_{3}=d_{2}\\ 3:-ru^{n+1}_{2}+(1+2r)u^{n+1}_{3}-ru^{n+1}_{4}=d_{3}\\ \cdots \\ i:-ru^{n+1}_{i-1}+(1+2r)u^{n+1}_{i}-ru^{n+1}_{i+1}=d_{i}\\ \cdots \\ N-2:-ru^{n+1}_{N-3}+(1+2r)u^{n+1}_{N-2}-ru^{n+1}_{N-1}=d_{N-2}\\ N-1:-ru^{n+1}_{N-2}+(1+2r)u^{n+1}_{N-1}=d_{N-1}+ru^{n+1}_{N}\\ \end{array}\end{split}\]

An illustration of a tridiagonal system is given in the following Equation.

\[\begin{split}\begin{bmatrix} b_{1}&c_{1} &0&\cdots &0 \\ a_{2}&b_{2}&c_{2} &\cdots &0 \\ \vdots&\vdots& & &\vdots \\ 0&\cdots&a_{N-2} &b_{N-2} &c_{N-2} \\ 0&0&\cdots &a_{N-1} &b_{N-1} \\ \end{bmatrix} \begin{bmatrix} u_{1}\\u_{2}\\\vdots \\u_{N-2}\\u_{N-1} \end{bmatrix} =\begin{bmatrix} \hat{d}_{1}\\ \hat{d}_{2}\\ \vdots \\ \hat{d}_{N-2}\\ \hat{d}_{N-1} \end{bmatrix}\end{split}\]

\[\begin{split}\begin{bmatrix} 1+2r&-r &0&\cdots &0 \\ -r&1+2r&-r &\cdots &0 \\ \vdots&\vdots& & &\vdots \\ 0&\cdots&-r &1+2r &-r \\ 0&0&\cdots &-r &1+2r \\ \end{bmatrix} \begin{bmatrix} u_{1}\\u_{2}\\\vdots \\u_{N-2}\\u_{N-1} \end{bmatrix} =\begin{bmatrix} \hat{d}_{1}\\ \hat{d}_{2}\\ \vdots \\ \hat{d}_{N-2}\\ \hat{d}_{N-1} \end{bmatrix}\end{split}\]

\[\begin{split}\begin{array}{l} a_{1}=0,b_{1}=1+2r,c_{1}=-r\\ a_{2}=-r,b_{2}=1+2r,c_{2}=-r\\ a_{3}=-r,b_{3}=1+2r,c_{3}=-r\\ \cdots\\ a_{i}=-r,b_{i}=1+2r,c_{i}=-r\\ \cdots\\ a_{N-2}=-r,b_{N-2}=1+2r,c_{N-2}=-r\\ a_{N-1}=-r,b_{N-1}=1+2r,c_{N-1}=0\\ \end{array}\end{split}\]

\[\begin{split}\begin{bmatrix} \hat{d}_{1}\\ \hat{d}_{2}\\ \vdots \\ \hat{d}_{N-2}\\ \hat{d}_{N-1} \end{bmatrix} =\begin{bmatrix} {d}_{1}+ru^{n+1}_{0}\\ {d}_{2}\\ \vdots \\ {d}_{N-2}\\ {d}_{N-1}+ru^{n+1}_{N} \end{bmatrix} =\begin{bmatrix} r{u}^{n}_{0}+(1-2r){u}^{n}_{1}+r{u}^{n}_{2}+ru^{n+1}_{0}\\ r{u}^{n}_{1}+(1-2r){u}^{n}_{2}+r{u}^{n}_{3}\\ \vdots \\ r{u}^{n}_{N-3}+(1-2r){u}^{n}_{N-2}+r{u}^{n}_{N-1}\\ r{u}^{n}_{N-2}+(1-2r){u}^{n}_{N-1}+r{u}^{n}_{N}+ru^{n+1}_{N} \end{bmatrix}\end{split}\]

\[\begin{split}\begin{bmatrix} {a}_{1}\\ {a}_{2}\\ \vdots \\ {a}_{N-2}\\ {a}_{N-1} \end{bmatrix} = \begin{bmatrix} {a}[0]\\ {a}[1]\\ \vdots \\ {a}[N-3]\\ {a}[N-2] \end{bmatrix} = \begin{bmatrix} 0\\ -r\\ \vdots \\ -r\\ -r \end{bmatrix}\end{split}\]

\[\begin{split}\begin{bmatrix} {b}_{1}\\ {b}_{2}\\ \vdots \\ {b}_{N-2}\\ {b}_{N-1} \end{bmatrix} = \begin{bmatrix} {b}[0]\\ {b}[1]\\ \vdots \\ {b}[N-3]\\ {b}[N-2] \end{bmatrix} = \begin{bmatrix} 1+2r\\ 1+2r\\ \vdots \\ 1+2r\\ 1+2r \end{bmatrix}\end{split}\]

\[\begin{split}\begin{bmatrix} {c}_{1}\\ {c}_{2}\\ \vdots \\ {c}_{N-2}\\ {c}_{N-1} \end{bmatrix} = \begin{bmatrix} {c}[0]\\ {c}[1]\\ \vdots \\ {c}[N-3]\\ {c}[N-2] \end{bmatrix} = \begin{bmatrix} -r\\ -r\\ \vdots \\ -r\\ 0 \end{bmatrix}\end{split}\]

The Thomas algorithm used for solving the tridiagonal matrix is ：

\[\begin{split}\begin{array}{l} 1:\text{Given a,b,c,d}\\ 2:\text{Allocate q}\quad\quad\diamondsuit \text{Storage of superdiagonal array }\\ 3:u_1=d_1/b_{1}\\ 4:\text{ for i=2 to N do}\quad\quad\diamondsuit\text{Forward elimination }\\ 5:\quad q_{i}=c_{i-1}/b_{i-1}\\ 6:\quad b_{i}=b_{i}-q_{i}a_{i}\\ 7:\quad u_{i}=(d_{i}-a_{i}u_{i-1})/b_{i}\\ 8:\text{end for}\\ 9:\text{ for i=N-1 to 1 do}\\ 10:\quad u_{i}=(u_{i}-q_{i+1}u_{i+1})\quad\quad\diamondsuit\text{Backward substitution }\\ 11:\text{end for}\\ \end{array}\end{split}\]

Matrix Configuration

Tri-diagonal systems for \(nx+1\) unknowns may be written as:

\[a_{i}u_{i-1}+b_{i}u_{i}+c_{i}u_{i+1}=d_{i}\]

Boundary points: \((1,nx+1)\), inner points: \((2,3,\cdots,nx-1,nx)\)

We know the values at the boundaries (\(B\)):

\[\begin{split}\begin{align} u_{1} & = B_{1}\\ u_{nx+1} & = B_{nx+1}\\ \end{align}\end{split}\]

the coefficient at the boundaries (\(B\)):

\[\begin{split}\begin{align} a_{1} & = 0\\ b_{1} & = 1\\ c_{1} & = 0\\ a_{nx+1} & = 0\\ b_{nx+1} & = 1\\ c_{nx+1} & = 0\\ \end{align}\end{split}\]

\[\begin{split}\begin{array}{l} a_{1}u_{0}+b_{1}u_{1}+c_{1}u_{2}&=0u_{0}+u_{1}+0u_{2} = u_{1}=B_{1}\\ a_{nx+1}u_{nx}+b_{nx+1}u_{nx+1}+c_{nx+1}u_{nx+2}&=0u_{nx}+u_{nx+1}+0u_{nx+2} = u_{nx+1}=B_{nx+1}\\ \end{array}\end{split}\]

So the matrix looks like this, with known coefficients \(a,b,c,d\). The vector \(u\) is unknown.

\[\begin{split}\left[\begin{array}{cccccc} 1 & 0 & 0 & \cdots & & 0 \\ a_{2} & b_{2} & c_{2} & & & \\ 0 & a_{3} & b_{3} & c_{3} & & \\ \vdots & & \ddots & \ddots & \ddots & \vdots \\ & & a_{n x-1} & b_{n x-1} & c_{n x-1} & 0 \\ 0 & & & a_{n x} & b_{n x} & c_{n x} \\ & & \cdots & 0 & 0 & 1 \end{array}\right]\left[\begin{array}{c} u_{1} \\ u_{2} \\ u_{3} \\ \vdots \\ u_{n x-1} \\ u_{n x} \\ u_{n x+1} \end{array}\right]=\left[\begin{array}{c} B_{1} \\ d_{2} \\ d_{3} \\ \vdots \\ d_{n x-1} \\ d_{n x} \\ B_{n x+1} \end{array}\right]\end{split}\]

Subscripts count from zero

\[\begin{split}\left[\begin{array}{cccccc} 1 & 0 & 0 & \cdots & & 0 \\ a_{1} & b_{1} & c_{1} & & & \\ 0 & a_{2} & b_{2} & c_{2} & & \\ \vdots & & \ddots & \ddots & \ddots & \vdots \\ & & a_{n x-2} & b_{n x-2} & c_{n x-2} & 0 \\ 0 & & & a_{n x-1} & b_{n x-1} & c_{n x-1} \\ & & \cdots & 0 & 0 & 1 \end{array}\right]\left[\begin{array}{c} u_{0} \\ u_{1} \\ u_{2} \\ \vdots \\ u_{n x-2} \\ u_{n x-1} \\ u_{n x} \end{array}\right]=\left[\begin{array}{c} B_{0} \\ d_{1} \\ d_{2} \\ \vdots \\ d_{n x-2} \\ d_{n x-1} \\ B_{n x} \end{array}\right]\end{split}\]

\[\begin{split}\left[\begin{array}{cccccc} b_{0} & c_{0} & 0 & \cdots & & 0 \\ a_{1} & b_{1} & c_{1} & & & \\ 0 & a_{2} & b_{2} & c_{2} & & \\ \vdots & & \ddots & \ddots & \ddots & \vdots \\ & & a_{n x-2} & b_{n x-2} & c_{n x-2} & 0 \\ 0 & & & a_{n x-1} & b_{n x-1} & c_{n x-1} \\ & & \cdots & 0 & a_{n x} & b_{n x} \end{array}\right]\left[\begin{array}{c} u_{0} \\ u_{1} \\ u_{2} \\ \vdots \\ u_{n x-2} \\ u_{n x-1} \\ u_{n x} \end{array}\right]=\left[\begin{array}{c} B_{0} \\ d_{1} \\ d_{2} \\ \vdots \\ d_{n x-2} \\ d_{n x-1} \\ B_{n x} \end{array}\right]\end{split}\]

\[a_{1} = 0,\quad b_{1} = 1, \quad c_{1}= 0\]

\[a_{nx+1} = 0,\quad b_{nx+1}=1,\quad c_{nx+1} = 0\]

1D heat equation: Forward time central space (FTCS) scheme
1D heat equation: Runge-Kutta Numerical Scheme