对于微分方程
\[ \frac{\partial Q}{\partial t} = R(Q) \]
BDF1和BDF2的理论公式如下
\[ \begin{aligned} Q^{n+1} - Q^{n} &= \Delta t R^{n+1} \\ Q^{n+1} - \frac{4}{3} Q^{n} + \frac{1}{3} Q^{n-1} &= \frac{2}{3} \Delta t R^{n+1} \\ \end{aligned} \]
BLUSGS算法描述见：Sun, Y.; Wang, Z. & Liu, Y. Efficient Implicit Non-linear LU-SGS Approach for Compressible Flow Computation Using High-order Spectral Difference Method, Comm. Comput. Phys., 2009, 5, 760-778
BDF1的原始形式
\[ \begin{aligned} \frac{ Q^{n+1} - Q^{n} }{\Delta t} &= R^{n+1} \\ \frac{ Q^{n+1} - Q^{n} }{\Delta t} - \left ( R^{n+1} - R^{n} \right ) &= R^{n} \\ \end{aligned} \]
线化\(R^{n+1}-R^{n}\)
\[ \begin{aligned} R^{n+1} - R^{n} = \frac{\partial R_c}{\partial Q_c} \Delta Q^{n+1}_c + \sum \frac{\partial R_{c}}{\partial Q_{nb}} \Delta Q^{n+1}_{nb} \end{aligned} \]
将线化\(R\)代入原始形式
\[ \begin{aligned} \Delta Q^{n+1} &= Q^{n+1} - Q^{n} \\ \frac{ \Delta Q^{n+1} }{\Delta t} - \frac{\partial R_c}{\partial Q_c} \Delta Q^{n+1}_c - \sum \frac{\partial R_{c}}{\partial Q_{nb}} \Delta Q^{n+1}_{nb} &= R^n_c \\ \left ( \frac{ I }{\Delta t} - \frac{\partial R_c}{\partial Q_c} \right ) \Delta Q^{n+1}_c - \sum \frac{\partial R_{c}}{\partial Q_{nb}} \Delta Q^{n+1}_{nb} &= R^n_c \\ \end{aligned} \]
Gauss-Seidel加速
\[ \begin{aligned} \left ( \frac{ I }{\Delta t} - \frac{\partial R_c}{\partial Q_c} \right ) \Delta Q^{k+1}_c - \sum \frac{\partial R_{c}}{\partial Q_{nb}} \Delta Q^{k}_{nb} &= R^n_c \\ \Delta Q^{k+1}_c &= Q^{k+1} - Q^{n} \\ \end{aligned} \]
再次线化\(R^{k}-R^{n}\)
\[ R^{k} - R^{n} = \frac{\partial R_c}{\partial Q_c} \Delta Q^{k}_c + \sum \frac{\partial R_{c}}{\partial Q_{nb}} \Delta Q^{k}_{nb} \]
将线化\(R\)代入，整理
\[ \begin{aligned} \left ( \frac{ I }{\Delta t} - \frac{\partial R_c}{\partial Q_c} \right ) \Delta Q^{k+1}_c &= R^k_c - \frac{\partial R_{c}}{\partial Q_{c}} \Delta Q^{k}_{c} \\ &= R^k_c + \left ( \frac{ I }{\Delta t} - \frac{\partial R_c}{\partial Q_c} \right ) \Delta Q^{k}_{c} - \frac{ I }{\Delta t} \Delta Q^{k}_{c} \\ \left ( \frac{ I }{\Delta t} - \frac{\partial R_c}{\partial Q_c} \right ) \left ( \Delta Q^{k+1}_c - \Delta Q^{k}_{c} \right ) &= R^k_c - \frac{ I }{\Delta t} \Delta Q^{k}_{c} \\ \left ( \frac{ I }{\Delta t} - \frac{\partial R_c}{\partial Q_c} \right ) \tilde{\Delta} Q^{k+1}_c &= R^k_c - \frac{ I }{\Delta t} \Delta Q^{k}_{c} \\ \tilde{\Delta} Q^{k+1}_c &= Q^{k+1}_c - Q^{k}_c \\ \end{aligned} \]
最后两个方程是最终的BDF1的BLUSGS实现。
BDF2的实现与BDF1类似。
BDF2的原始形式
\[ \begin{aligned} \frac{ Q^{n+1} - \frac{4}{3} Q^{n} + \frac{1}{3} Q^{n-1} }{ \frac{2}{3} \Delta t } &= R^{n+1} \\ \frac{3}{2} \left ( \frac{ \Delta Q^{n+1} }{ \Delta t } - \frac{1}{3} \frac{ \Delta Q^{n} }{ \Delta t } \right ) - \left ( R^{n+1} - R^{n} \right ) &= R^{n} \\ \Delta Q^{n+1} &= Q^{n+1} - Q^{n} \\ \end{aligned} \]
线化\(R^{n+1}-R^{n}\)
\[ \begin{aligned} R^{n+1} - R^{n} = \frac{\partial R_c}{\partial Q_c} \Delta Q^{n+1}_c + \sum \frac{\partial R_{c}}{\partial Q_{nb}} \Delta Q^{n+1}_{nb} \end{aligned} \]
将线化\(R\)代入原始形式
\[ \begin{aligned} \left ( \frac{3}{2} \frac{ I }{\Delta t} - \frac{\partial R_c}{\partial Q_c} \right ) \Delta Q^{n+1}_c - \frac{1}{2} \frac{ I }{\Delta t} \Delta Q^{n}_c - \sum \frac{\partial R_{c}}{\partial Q_{nb}} \Delta Q^{n+1}_{nb} &= R^n_c \\ \end{aligned} \]
Gauss-Seidel加速
\[ \begin{aligned} \left ( \frac{3}{2} \frac{ I }{\Delta t} - \frac{\partial R_c}{\partial Q_c} \right ) \Delta Q^{k+1}_c - \frac{1}{2} \frac{ I }{\Delta t} \Delta Q^{n}_c - \sum \frac{\partial R_{c}}{\partial Q_{nb}} \Delta Q^{k}_{nb} &= R^n_c \\ \Delta Q^{k+1}_c &= Q^{k+1} - Q^{n} \\ \end{aligned} \]
再次线化\(R^{k}-R^{n}\)
\[ R^{k} - R^{n} = \frac{\partial R_c}{\partial Q_c} \Delta Q^{k}_c + \sum \frac{\partial R_{c}}{\partial Q_{nb}} \Delta Q^{k}_{nb} \]
将线化\(R\)代入，整理
\[ \begin{aligned} \left ( \frac{3}{2} \frac{ I }{\Delta t} - \frac{\partial R_c}{\partial Q_c} \right ) \Delta Q^{k+1}_c - \frac{1}{2} \frac{ I }{\Delta t} \Delta Q^{n}_c &= R^k_c - \frac{\partial R_{c}}{\partial Q_{c}} \Delta Q^{k}_{c} \\ &= R^k_c + \left ( \frac{3}{2} \frac{ I }{\Delta t} - \frac{\partial R_c}{\partial Q_c} \right ) \Delta Q^{k}_{c} - \frac{3}{2} \frac{ I }{\Delta t} \Delta Q^{k}_{c} \\ \left ( \frac{3}{2} \frac{ I }{\Delta t} - \frac{\partial R_c}{\partial Q_c} \right ) \left ( \Delta Q^{k+1}_c - \Delta Q^{k}_{c} \right ) - \frac{1}{2} \frac{ I }{\Delta t} \Delta Q^{n}_c &= R^k_c - \frac{3}{2} \frac{ I }{\Delta t} \Delta Q^{k}_{c} \\ \left ( \frac{3}{2} \frac{ I }{\Delta t} - \frac{\partial R_c}{\partial Q_c} \right ) \tilde{\Delta} Q^{k+1}_c &= R^k_c - \frac{3}{2} \frac{ I }{\Delta t} \Delta Q^{k}_{c} + \frac{1}{2} \frac{ I }{\Delta t} \Delta Q^{n}_c \\ \left ( \frac{3}{2} \frac{ I }{\Delta t} - \frac{\partial R_c}{\partial Q_c} \right ) \tilde{\Delta} Q^{k+1}_c &= R^k_c - \frac{1}{2 \Delta t} \left ( 3 I Q^{k}_{c} - 4 I Q^{n}_c + I Q^{n-1}_c \right ) \\ \tilde{\Delta} Q^{k+1}_c &= Q^{k+1}_c - Q^{k}_c \\ \end{aligned} \]
最后两个方程是最终的BDF2的BLUSGS实现。
The BDF2OPT schemes are referred to \(BDF2OPT(\beta)\) which is defined in the following paper,
Vatsa, V., Carpenter, M., & Lockard, D. (2010). Re-evaluation of an optimized second-order backward difference (BDF2OPT) scheme for unsteady flow applications. In 48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition (p. 122).
Basically, BDF2OPT is a combination of BDF2 and BDF3 or higher-order BDF schemes, i.e.
\[ \text{BDF2OPT} = \beta \text{BDF3} + (1-\beta) \text{BDF2} \]
To ensure it to be A-stable, \(\beta \leq 0.5\). The paper recommends \(\beta=0.48\). In REF, BDF2OPT takes \(\beta=0.55\) and \([A,B,C] = [10/6, 5/6, 1/6]\).
Its BLUSGS implementation is the following procedure.
First, rewrite the combination to be as follows,
\[ \frac{\Delta Q^{n+1} - \frac{1}{3} \Delta Q^{n} - \frac{10}{33} \beta \Delta Q^{n} + \frac{6}{33} \beta \Delta Q^{n-1}}{ \left ( \frac{2}{3} - \frac{4}{33} \beta \right ) \Delta t } - \left ( R^{n+1} - R^{n} \right ) = R^{n} \]
where \(\Delta Q^{n+1} = Q^{n+1} - Q^{n}\).
Split the equation in terms of \(\Delta Q\) as follows,
\[ \begin{aligned} \frac{A}{\Delta t} \Delta Q^{n+1} - \frac{B}{\Delta t} \Delta Q^{n} &+ \frac{C}{\Delta t} \Delta Q^{n-1} - \left ( R^{n+1} - R^{n} \right ) = R^{n} \\ A &= \frac{1}{\frac{2}{3} - \frac{4}{33} \beta} \\ B &= \frac{\frac{1}{3} + \frac{10}{33} \beta}{\frac{2}{3} - \frac{4}{33} \beta} \\ C &= \frac{\frac{6}{33} \beta}{\frac{2}{3} - \frac{4}{33} \beta} \\ \end{aligned} \]
Substitute the linearized \(R\) which is
\[ R^{n+1} - R^{n} = \frac{\partial R_c}{\partial Q_c} \Delta Q^{n+1}_c + \sum \frac{\partial R_{c}}{\partial Q_{nb}} \Delta Q^{n+1}_{nb} \]
to obtain
\[ \left ( \frac{A}{\Delta t} I - \frac{\partial R_c}{\partial Q_c} \right ) \Delta Q^{n+1}_c - \frac{B}{\Delta t} \Delta Q^{n}_c + \frac{C}{\Delta t} \Delta Q^{n-1}_c - \sum \frac{\partial R_{c}}{\partial Q_{nb}} \Delta Q^{n+1}_{nb} = R^{n}_c \]
Employ Gauss-Seidel acceleration to obtain
\[ \left ( \frac{A}{\Delta t} I - \frac{\partial R_c}{\partial Q_c} \right ) \Delta Q^{k+1}_c - \frac{B}{\Delta t} \Delta Q^{n}_c + \frac{C}{\Delta t} \Delta Q^{n-1}_c - \sum \frac{\partial R_{c}}{\partial Q_{nb}} \Delta Q^{k}_{nb} = R^{n}_c \]
where \(\Delta Q^{k+1} = Q^{k+1} - Q^{n}\).
Substitute the linearized \(R^{k}_c - R^{n}_c\),
\[ R^{k}_c - R^{n}_c = \frac{\partial R_c}{\partial Q_c} \Delta Q^{k}_c + \sum \frac{\partial R_{c}}{\partial Q_{nb}} \Delta Q^{k}_{nb} \]
to obtain
\[ \left ( \frac{A}{\Delta t} I - \frac{\partial R_c}{\partial Q_c} \right ) \Delta Q^{k+1}_c - \frac{B}{\Delta t} \Delta Q^{n}_c + \frac{C}{\Delta t} \Delta Q^{n-1}_c = R^k_c - \frac{\partial R_{c}}{\partial Q_{c}} \Delta Q^{k}_{c} \]
Manipulate the equation to be as follows,
\[ \begin{aligned} \left ( \frac{A}{\Delta t} I - \frac{\partial R_c}{\partial Q_c} \right ) \tilde{\Delta} Q^{k+1}_c &= R^k_c - \frac{A}{\Delta t} \Delta Q^{k}_c + \frac{B}{\Delta t} \Delta Q^{n}_c - \frac{C}{\Delta t} \Delta Q^{n-1}_c \\ A &= \frac{1}{\frac{2}{3} - \frac{4}{33} \beta} \\ B &= \frac{\frac{1}{3} + \frac{10}{33} \beta}{\frac{2}{3} - \frac{4}{33} \beta} \\ C &= \frac{\frac{6}{33} \beta}{\frac{2}{3} - \frac{4}{33} \beta} \\ \end{aligned} \]
where \(\tilde{\Delta} Q^{k+1} = Q^{k+1} - Q^{k}\).
Now we have the BLUSGS implementation of the BDF2OPT.
Following this notation, BDF2 has the following form
\[ \begin{aligned} \left ( \frac{A}{\Delta t} I - \frac{\partial R_c}{\partial Q_c} \right ) \tilde{\Delta} Q^{k+1}_c &= R^k_c - \frac{A}{\Delta t} \Delta Q^{k}_c + \frac{B}{\Delta t} \Delta Q^{n}_c \\ A &= 1.5 \\ B &= 0.5 \\ \end{aligned} \]
BDF1 has the following form
\[ \begin{aligned} \left ( \frac{A}{\Delta t} I - \frac{\partial R_c}{\partial Q_c} \right ) \tilde{\Delta} Q^{k+1}_c &= R^k_c - \frac{A}{\Delta t} \Delta Q^{k}_c \\ A &= 1 \\ \end{aligned} \]
Refer to Holst’s JCP paper, Kennedy’s review, and Persson’s paper.
Christopher A. Kennedy, and Mark H. Carpenter. Diagonally implicit Runge-Kutta methods for ordinary differential equations. A review. NASA/TM-2016-219173(2016).
K.R. Holst et al., On the effect of temporal error in high-order simulations of unsteady flows, J. Comput. Phys. (2019), 108989, DOI: https://doi.org/10.1016/j.jcp.2019.108989.
Persson Per-Olof, High-order LES simulations using implicit-explicit Runge-Kutta schemes, 49th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition(2011), DOI: https://doi.org/10.2514/6.2011-684.
RK methods has a unified formula as follows,
\[ \begin{aligned} y_{n+1} &= y_{n}+h \sum_{i=1}^{s} b_{i} k_{i}, \\ k_{i} &= f \left( t_{n}+c_{i} h, y_{n}+h \sum_{j=1}^{s} a_{i j} k_{j} \right), \quad i=1, \ldots, s \\ \end{aligned} \]
The coefficients are stored in the Butcher tableau.
\[ \begin{array}{c|cccc} c_{1} & a_{11} & a_{12} & \ldots & a_{1 s} \\ c_{2} & a_{21} & a_{22} & \ldots & a_{2 s} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ c_{s} & a_{s 1} & a_{s 2} & \ldots & a_{s s} \\ \hline & b_{1} & b_{2} & \ldots & b_{s} \end{array} = \begin{array}{c|c} \mathbf{c} & \mathbf{A} \\ \hline & \mathbf{b^{T}} \\ \end{array} \]
where \(A\) is the RK matrix, \(b\) is the weights vector and \(c\) is the nodes vector.
RK methods may be either explicit (ERK) or implicit (IRK), and a particular method can be categorized using the RK matrix. If the matrix is only lower triangular, with the main diagonal zero, the method is explicit. If any value along the main diagonal or above it is non-zero, the method is implicit. Explicit schemes will not be discussed here because they are not A-stable.
There are several categories within the family of IRK methods.
SDIRK33 is taken from Persson’s paper. The first 3 in
the name is 3rd order and the second 3 in the name is 3
stages. The Butcher tabuleau is as follows,
\[ \begin{array}{c|c} \mathbf{c} & \mathbf{A} \\ \hline & \mathbf{b^{T}} \\ \end{array} = \begin{array}{c|ccc} 0.4358665215 & 0.4358665215 & 0 & 0 \\ 0.7179332608 & 0.2820667392 & 0.4358665215 & 0 \\ 1 & 1.208496649 & -0.644363171 & 0.4358665215 \\ \hline & 1.208496649 & -0.644363171 & 0.4358665215 \end{array} \]
The above matrix coefficients are in fact originally as follows,
\[ \begin{array}{c|ccc} \gamma & \gamma & 0 & 0 \\ c_{2} & \left(c_{2}-\gamma\right) & \gamma & 0 \\ 1 & \left(1-b_{2}-\gamma\right) & b_{2} & \gamma \\ \hline & \left(1-b_{2}-\gamma\right) & b_{2} & \gamma \end{array} \]
where the coefficients are as follows,
\[ \begin{aligned} b_{2} &= \frac{-3 \alpha^{2}}{4 \beta} \\ c_{2} &= \frac{2-9 \gamma+6 \gamma^{2}}{3 \alpha} \\ \alpha &= \left(1-4 \gamma+2 \gamma^{2}\right) \\ \beta &= \left(-1+6 \gamma-9 \gamma^{2}+3 \gamma^{3}\right) \\ \end{aligned} \]
\(\gamma = 0.43586652150845899941601945\) gives the L-stable scheme. Here is a Python script to show the calculation as follows,
from mpmath import mp, mpf, matrix
mp.dps = 32
gamma = mpf('0.43586652150845899941601945')
alpha = 1 - 4*gamma + 2*gamma**2
beta = -1 + 6*gamma - 9*gamma**2 + 3*gamma**3
b2 = (-3*alpha**2) / (4*beta)
c2 = (2-9*gamma+6*gamma**2) / (3*alpha)
A = matrix(3,3)
A[0,0] = gamma
A[1,1] = gamma
A[2,2] = gamma
A[1,0] = c2 - gamma
A[2,0] = 1 - b2 - gamma
A[2,1] = b2
print("A is ")
print(A)
B = A[2,:]
print("B is ")
print(B)
C = matrix(3,1)
C[0,0] = gamma
C[1,0] = c2
C[2,0] = 1
print("C is ")
print(C)The high precision version of the Butcher tableau is as follows,
A is
[ 0.43586652150845899941601945 0.0 0.0]
[0.28206673924577050029199027679033 0.43586652150845899941601945 0.0]
[ 1.2084966491760100703364776750294 -0.64436317068446906975249712502944 0.43586652150845899941601945]
B is
[1.2084966491760100703364776750294 -0.64436317068446906975249712502944 0.43586652150845899941601945]
C is
[ 0.43586652150845899941601945]
[0.71793326075422949970800972679033]
[ 1.0]Now we have the formula. How could we solve it? At the first stage of the SDIRK33, the formula is actually as follows,
\[ k_{1} = f \left( t_{n}+c_{1} h, y_{n} + h a_{11} k_{1} \right) \]
This is a nonlinear equation. We could manipulate it to adapt to the LUSGS solver.
\[ \begin{aligned} \frac{ y_{n1} - y_{n} }{ h a_{11} } & = k_{1} \\ k_{1} &= f \left( t_{n1}, y_{n1} \right) \\ y_{n1} &= y_{n} + h a_{11} k_{1} \\ t_{n1} &= t_{n}+c_{1} h \\ \end{aligned} \]
This is very similar to BDF1. Probably this is why the stage order is two.
Similarly, at the stage no. 2, the formula is manipulated to adapt to the LUSGS solver as follows,
\[ \begin{aligned} k_{2} &= f \left( t_{n}+c_{2} h, y_{n} + h a_{21} k_{1} + h a_{22} k_{2} \right) \\ \frac{ y_{n2} - \left ( y_{n} + h a_{21} k_{1} \right ) }{ h a_{22} } &= k_{2} \\ k_{2} &= f \left( t_{n2}, y_{n2} \right) \\ y_{n2} &= y_{n} + h a_{21} k_{1} + h a_{22} k_{2} \\ t_{n2} &= t_{n}+c_{2} h \\ \end{aligned} \]
And at the stage no. 3, the formula is manipulated to adapt to the LUSGS solver as follows,
\[ \begin{aligned} k_{3} &= f \left( t_{n}+c_{3} h, y_{n} + h a_{31} k_{1} + h a_{32} k_{2} + h a_{33} k_{3} \right) \\ \frac{ y_{n3} - \left ( y_{n} + h a_{31} k_{1} + h a_{32} k_{2} \right ) }{ h a_{33} } &= k_{3} \\ k_{3} &= f \left( t_{n3}, y_{n3} \right) \\ y_{n3} &= y_{n} + h a_{31} k_{1} + h a_{32} k_{2} + h a_{33} k_{3} \\ t_{n3} &= t_{n}+c_{3} h \\ \end{aligned} \]