Time-Marching Schemes

Problem Setup

We focus on integrating several partial differential equations involving both temporal and spatial derivatives in this project. When focusing only on the temporal derivatives, all equations can be written as

\[\newcommand{\oldvar}[1]{{#1^{n}}} \newcommand{\newvar}[1]{{#1^{n+1}}} \newcommand{\oerror}[1]{\mathcal{O} \left( \Delta t^{#1} \right)} \frac{df}{dt} = g.\]

Note

We assume the right-hand-side term \(g\) is not a function of time explicitly, i.e., the following discussion is limited for autonomous systems. This applies to all equations discussed in this project.

This is equivalent to

\[\int_\oldvar{t}^\newvar{t} \frac{df}{dt} dt = \int_\oldvar{t}^\newvar{t} g dt,\]

or

\[\newvar{f} - \oldvar{f} = \int_\oldvar{t}^\newvar{t} g dt,\]

where we introduce

\[ \begin{align}\begin{aligned}\oldvar{f} & \equiv f \left( t = \oldvar{t} \right),\\\newvar{f} & \equiv f \left( t = \newvar{t} \right),\end{aligned}\end{align} \]

and

\[\Delta t \equiv \newvar{t} - \oldvar{t}.\]

To obtain \(\newvar{f}\) from \(\oldvar{f}\), we need to evaluate the right-hand-side integral. In this page, we aim at approximating the right-hand-side term using the Taylor-series expansion around \(\oldvar{t}\):

\[ \begin{align}\begin{aligned}\newvar{f} & = \oldvar{f} + \frac{d\oldvar{f}}{dt} \Delta t + \frac{1}{2} \frac{d^2\oldvar{f}}{dt^2} \left( \Delta t \right)^2 + \oerror{3}\\& = \oldvar{f} + \oldvar{g} \Delta t + \frac{1}{2} \frac{d\oldvar{g}}{dt} \left( \Delta t \right)^2 + \oerror{3}.\end{aligned}\end{align} \]

Zeroth-Order Scheme

The simplest but useless way to approximate the above Taylor series expansion is

\[\newvar{f} \approx \oldvar{f},\]

where we only pick-up the first term. The leading-order error (local truncation error) is

\[\oldvar{g} \Delta t = \oerror{1},\]

which indicates the first-order accuracy. Since our objective is integrate the equation for a long time, we should consider the global truncation error, giving a reduced order of accuracy: \(\oerror{0}\). This indicates that the error never shrinks even with a smaller \(\Delta t\), which is useless.

Note

Hereafter we only focus on the global truncation errors.

Euler-Forward Scheme

We obtain a simple but (to some extent) practical scheme by extracting the first two terms:

\[\newvar{f} = \oldvar{f} + \oldvar{g} \Delta t,\]

or equivalently

\[ \begin{align}\begin{aligned}\Delta f & = \oldvar{g} \Delta t,\\\newvar{f} & = \oldvar{f} + \Delta f,\end{aligned}\end{align} \]

which is known as the Euler-forward (or Euler-explicit) scheme.

Since the deviation from the series expansion is \(\oerror{2}\) (locally), this scheme has the first-order accuracy (globally).

Although it is very simple and easy to use, two clear issues exist:

• It only has the first-order accuracy in time,

• All terms are treated explicitly.

Fully-Explicit Runge-Kutta Scheme

To achieve a higher-order accuracy, we need to extract more than two terms from the series expansion above. In this project, we adopt the explicit Runge-Kutta scheme. Since the system is autonomous, we use a three-step and low-storage scheme to achieve the second-order accuracy in time:

\[ \begin{align}\begin{aligned}& \text{do}\,\,k = 0, 2\\& \,\,\,\, \Delta f = \left( \alpha^k g^{k } + \beta^k g^{k-1} \right) \Delta t,\\& \,\,\,\, f^{k+1} = f^{k } + \Delta f.\\& \text{enddo}\end{aligned}\end{align} \]

Although the coefficients \(\alpha^k,\beta^k\) have several possibilities, we adopt

\[\left( \alpha^0, \alpha^1, \alpha^2 \right) = \left( \frac{32}{60}, \frac{25}{60}, \frac{45}{60} \right),\]

and

\[\left( \beta^0, \beta^1, \beta^2 \right) = \left( \frac{0}{60}, -\frac{17}{60}, -\frac{25}{60} \right),\]

following e.g., Rai and Moin, J. Comput. Phys. (96), 1991, Verzicco and Orlandi, J. Comput. Phys. (123), 1996, Costa, Comput. Math Appl. (76), 2018.

For later convenience, we also introduce

\[\gamma^{k } \equiv \alpha^{k } + \beta^{k },\]

or explicitly

\[\left( \gamma^0, \gamma^1, \gamma^2 \right) = \left( \frac{32}{60}, \frac{8}{60}, \frac{20}{60} \right).\]

Note that, since we have

\[\sum_{k = 0}^{2} \gamma^k = 1,\]

the above three-step Runge-Kutta scheme is essentially a combination of three Euler-forward scheme with \(\gamma^k \Delta t\) as time-step sizes.

src/runge_kutta.c

static const double a0 = + 32. / 60., b0 =    0. / 60.;
static const double a1 = + 25. / 60., b1 = - 17. / 60.;
static const double a2 = + 45. / 60., b2 = - 25. / 60.;
const rkcoef_t rkcoefs[RKSTEPMAX] = {
  {.alpha = a0, .beta = b0, .gamma = a0 + b0},
  {.alpha = a1, .beta = b1, .gamma = a1 + b1},
  {.alpha = a2, .beta = b2, .gamma = a2 + b2},
};

Proof of the second-order accuracy

Because of

\[f^{k+1} = f^{k } + \alpha^k g^{k } \Delta t + \beta^k g^{k-1} \Delta t,\]

we find

\[g^{k+1} = \frac{df^{k+1}}{dt} = \frac{df^k}{dt} + \alpha^k \frac{d^2f^k}{dt^2} \Delta t + \beta^k \frac{d^2f^{k-1}}{dt^2} \Delta t.\]

We use this relation repeatedly. For \(k = 1\), we have

\[f^1 = f^n + \alpha^0 g^n \Delta t = f^n + \alpha^0 \frac{df^n}{dt} \Delta t,\]

whose derivation leads to

\[g^1 = \frac{df^n}{dt} + \alpha^0 \frac{d^2f^n}{dt^2} \Delta t.\]

Note that \(k = 0\) corresponds to \(n\) (old information). For \(k = 2\), we have

\[ \begin{align}\begin{aligned}f^2 & = f^1 + \alpha^1 g^1 \Delta t + \beta^1 g^0 \Delta t\\& = f^n + \alpha^0 \frac{df^n}{dt} \Delta t + \alpha^1 \left( \frac{df^n}{dt} + \alpha^0 \frac{d^2f^n}{dt^2} \Delta t \right) \Delta t + \beta^1 \frac{df^n}{dt} \Delta t,\\& = f^n + \left( \alpha^0 + \alpha^1 + \beta^1 \right) \frac{df^n}{dt} \Delta t + \alpha^0 \alpha^1 \frac{d^2f^n}{dt^2} \Delta t^2,\end{aligned}\end{align} \]

whose derivation leads to

\[g^2 = \frac{df^n}{dt} + \left( \alpha^0 + \alpha^1 + \beta^1 \right) \frac{d^2f^n}{dt^2} \Delta t + \alpha^0 \alpha^1 \frac{d^3f^n}{dt^3} \Delta t^2.\]

For \(k = 3\), we have

\[ \begin{align}\begin{aligned}f^3 & = f^2 + \alpha^2 g^2 \Delta t + \beta^2 g^1 \Delta t\\& = f^n + \left( \alpha^0 + \alpha^1 + \beta^1 \right) \frac{df^n}{dt} \Delta t + \alpha^0 \alpha^1 \frac{d^2f^n}{dt^2} \Delta t^2 + \alpha^2 \left( \frac{df^n}{dt} + \left( \alpha^0 + \alpha^1 + \beta^1 \right) \frac{d^2f^n}{dt^2} \Delta t + \alpha^0 \alpha^1 \frac{d^3f^n}{dt^3} \Delta t^2 \right) \Delta t + \beta^2 \left( \frac{df^n}{dt} + \alpha^0 \frac{d^2f^n}{dt^2} \Delta t \right) \Delta t,\\& = f^n + \left( \alpha^0 + \alpha^1 + \beta^1 + \alpha^2 + \beta^2 \right) \frac{df^n}{dt} \Delta t + \left( \alpha^0 \alpha^1 + \alpha^0 \alpha^2 + \alpha^1 \alpha^2 + \alpha^2 \beta^1 + \alpha^0 \beta^2 \right) \frac{d^2f^n}{dt^2} \Delta t^2 + \oerror{3}\\& = f^n + \frac{df^n}{dt} \Delta t + \frac{1}{2} \frac{d^2f^n}{dt^2} \Delta t^2 + \oerror{3}.\end{aligned}\end{align} \]

Since \(k = 3\) corresponds to \(n + 1\) (new information), we find

\[\newvar{f} = f^n + \frac{df^n}{dt} \Delta t + \frac{1}{2} \frac{d^2f^n}{dt^2} \left( \Delta t^2 \right) + \oerror{3},\]

indicating that this scheme has the third-order accuracy (locally) and the second-order accuracy (globally).

Implicit Treatment

As justified later, we treat some terms implicitly in time; among others we adopt Crank-Nicolson scheme:

\[ \begin{align}\begin{aligned}\Delta f & = \frac{1}{2} \oldvar{g} \Delta t + \frac{1}{2} g^{n+1} \Delta t,\\\newvar{f} - \oldvar{f} & = \Delta f,\end{aligned}\end{align} \]

which has the second-order accuracy in time as well.

Proof of the second-order accuracy

Because of

\[g^{n+1} = \frac{d\newvar{f}}{dt} = \frac{df^n}{dt} + \frac{d^2f^n}{dt^2} \Delta t + \frac{1}{2} \frac{d^3f^n}{dt^3} \Delta t^2 + \oerror{3},\]

we find

\[ \begin{align}\begin{aligned}\newvar{f} & = f^n + \frac{1}{2} \frac{df^n}{dt} \Delta t + \frac{1}{2} \left( \frac{df^n}{dt} + \frac{d^2f^n}{dt^2} \Delta t + \frac{1}{2} \frac{d^3f^n}{dt^3} \Delta t^2 + \frac{df^n}{dt} \right) \Delta t\\& = f^n + \frac{df^n}{dt} \Delta t + \frac{1}{2} \frac{d^2f^n}{dt^2} \Delta t^2 + \oerror{3},\end{aligned}\end{align} \]

indicating that this scheme has the third-order accuracy (locally) and the second-order accuracy (globally).

Combined scheme

We embed the implicit treatment (Crank-Nicolson scheme) into the Runge-Kutta iterations:

\[\begin{split}&\text{do}\,\, k = 0, 2 \\ &\,\,\,\, \Delta f = \left( \alpha^k g^{k } + \beta^k g^{k-1} \right) \Delta t + \frac{1}{2} \left( h^{k } + h^{k+1} \right) \gamma^k \Delta t, \\ &\,\,\,\, f^{k+1} = f^{k } + \Delta f, \\ &\text{enddo}\end{split}\]

where \(g\) and \(h\) are terms treated explicitly and implicitly in time, respectively.