ylee

The recursive system-free method

Posted on: Wed 16 October 2024

The system-free (SF) method [2] (see this article) demonstrated its efficiency and high-order accuracy by combining it with the Picard integration formulation (PIF) method. [3] By using the SF approach, the analytical derivations of flux Jacobians and Hessians can be approximated numerically without affecting the accuracy and stability of the PIF method. However, extending the third-order PIF method to a higher-order accuracy involves calculating higher degrees of Jacobian-like tensor-vector contractions (e.g., ∂3F∂U3 \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \pddd{\bF}{\bU}), significantly adding to the complexity of such schemes.

The original system-free approach

Following the same algorithmic strategy, the fourth-order extension of the PIF method requires the third-order derivatives of the flux function. Theoretically speaking, the original SF method can approximate this term as,

∂3F∂U3⋅V⋅V⋅V=12ε3[−F(U−2εV)+2F(U−εV)−2F(U+εV)+F(U+2εV)]+O(ε2), \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \begin{aligned} \pddd{\bF}{\bU} \cdot \bV \cdot \bV \cdot \bV = \frac{1}{2 \varepsilon^{3}} \bigg[ -\bF(\bU - 2 \varepsilon \bV) + 2\bF(\bU - \varepsilon \bV) - 2\bF(\bU + \varepsilon \bV)+ \mathbf{F}(\bU + 2 \varepsilon \bV) \bigg] + \mathcal{O} ( \varepsilon^{2} ), \end{aligned}

for the tensor contractions between the same vector, V \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \bV. In order to approximate the case of different vectors, (V \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \bV, W \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \bW, and X \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \bX for example), the above equation should be computed multiple times:

∂3F∂U3⋅V⋅W⋅X=16[∂3F∂U3⋅(V+W+X)⋅(V+W+X)⋅(V+W+X)−∂3F∂U3⋅(V+W)⋅(V+W)⋅(V+W)−∂3F∂U3⋅(V+X)⋅(V+X)⋅(V+X)−∂3F∂U3⋅(W+X)⋅(W+X)⋅(W+X)+∂3F∂U3⋅V⋅V⋅V+∂3F∂U3⋅W⋅W⋅W+∂3F∂U3⋅X⋅X⋅X]. \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \begin{aligned} \pddd{\bF}{\bU} \cdot \bV \cdot \bW \cdot \bX = \frac{1}{6} \bigg[ &\pddd{\bF}{\bU} \cdot \left( \bV + \bW + \bX \right) \cdot \left( \bV + \bW + \bX \right) \cdot \left( \bV + \bW + \bX \right) \\ & -\pddd{\bF}{\bU} \cdot \left( \bV + \bW \right) \cdot \left( \bV + \bW \right) \cdot \left( \bV + \bW \right) \\ & -\pddd{\bF}{\bU} \cdot \left( \bV + \bX \right) \cdot \left( \bV + \bX \right) \cdot \left( \bV + \bX \right) \\ & -\pddd{\bF}{\bU} \cdot \left( \bW + \bX \right) \cdot \left( \bW + \bX \right) \cdot \left( \bW + \bX \right) \\ & + \pddd{\bF}{\bU} \cdot \bV \cdot \bV \cdot \bV + \pddd{\bF}{\bU} \cdot \bW \cdot \bW \cdot \bW + \pddd{\bF}{\bU} \cdot \bX \cdot \bX \cdot \bX \bigg]. \end{aligned}

As shown above, the original SF method requires to repeat the ∂3F∂U3⋅V⋅V⋅V \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \pddd{\bF}{\bU} \cdot \bV \cdot \bV \cdot \bV approximation seven times just for getting a single term, ∂3F∂U3⋅V⋅W⋅X \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \pddd{\bF}{\bU} \cdot \bV \cdot \bW \cdot \bX. Although it is possible to use for the fourth-order PIF method, it is not very practical because of the code complexity and slow performance.

The recursive approach

The newly improved version of SF method is taking a "recursive" strategy for approximating high degree of Jacobian-like tensor terms. In order to take a recursive strategy, let's define a functional Du \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \mathcal{D}_{u} as,

∂F∂U⋅V≈Du(F ; V)≔12εv[F(U+εvV)−F(U−εvV)], \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \pd{\bF}{\bU} \cdot \bV \approx \mathcal{D}_{u} (\bF \scolon \bV) \coloneqq \frac{1}{2\varepsilon_{v}} \bigg[ \bF(\bU + \varepsilon_{v} \bV) - \bF(\bU - \varepsilon_{v} \bV) \bigg],

which is the Jacobian-free approximation from the original SF method. For higher-order derivatives, we take the functional Du \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \mathcal{D}_{u} in a successive fashion. For example, the Hessian-vector contraction would be,

∂2F∂U2⋅V⋅W≈Du(Du(F ; V) ; W)=14εvεw[F(U+εvV+εwW)−F(U−εvV+εwW)−F(U+εvV−εwW)+F(U−εvV−εwW)]. \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \begin{aligned} \pdd{\bF}{\bU} \cdot \bV \cdot \bW &\approx \mathcal{D}_{u} \Big( \mathcal{D}_{u} (\bF \scolon \bV) \scolon \bW \Big) \\ =\frac{1}{4 \varepsilon_{v} \varepsilon_{w}} \bigg[ &\bF(\bU + \varepsilon_{v} \bV + \varepsilon_{w} \bW) -\bF(\bU - \varepsilon_{v} \bV + \varepsilon_{w} \bW)\\ -&\bF(\bU + \varepsilon_{v} \bV - \varepsilon_{w} \bW) +\bF(\bU - \varepsilon_{v} \bV - \varepsilon_{w} \bW) \bigg]. \end{aligned}

Note that the above improved version of Hessian-free method is applicable regardless the tensor contraction with two identical vectors or with two distinct vectors, hence it does not require separate formulations needed in the original SF approach.

The simplicity gain from the improved version of the SF method is further rewarded when considering the higher-order derivatives of F \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \bF. Following the equivalent strategy, the tensor contraction of the third-order derivative of the flux function, ∂3F∂U3 \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \pddd{\bF}{\bU} with three distinct vectors, V,W \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \bV, \bW and X \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \bX can be expressed in a compact form as,

∂3F∂U3⋅V⋅W⋅X≈Du(Du(Du(F ; V) ; W) ; X)=18εvεwεx[F(U+εvV+εwW+εxX)−F(U−εvV+εwW+εxX)−F(U+εvV−εwW+εxX)+F(U−εvV−εwW+εxX)−F(U+εvV+εwW−εxX)+F(U−εvV+εwW−εxX)+F(U+εvV−εwW−εxX)−F(U−εvV−εwW−εxX)]. \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \begin{aligned} \pddd{\bF}{\bU} \cdot \bV \cdot \bW \cdot \bX &\approx \mathcal{D}_{u} \Bigg( \mathcal{D}_{u} \Big( \mathcal{D}_{u} (\bF \scolon \bV) \scolon \bW \Big) \scolon \bX \Bigg) \\ =\frac{1}{8 \varepsilon_{v} \varepsilon_{w} \varepsilon_{x}} \bigg[ &\bF(\bU + \varepsilon_{v} \bV + \varepsilon_{w} \bW + \varepsilon_{x} \bX) -\bF(\bU - \varepsilon_{v} \bV + \varepsilon_{w} \bW + \varepsilon_{x} \bX)\\ -&\bF(\bU + \varepsilon_{v} \bV - \varepsilon_{w} \bW + \varepsilon_{x} \bX) +\bF(\bU - \varepsilon_{v} \bV - \varepsilon_{w} \bW + \varepsilon_{x} \bX)\\ -&\bF(\bU + \varepsilon_{v} \bV + \varepsilon_{w} \bW - \varepsilon_{x} \bX) +\bF(\bU - \varepsilon_{v} \bV + \varepsilon_{w} \bW - \varepsilon_{x} \bX)\\ +&\bF(\bU + \varepsilon_{v} \bV - \varepsilon_{w} \bW - \varepsilon_{x} \bX) -\bF(\bU - \varepsilon_{v} \bV - \varepsilon_{w} \bW - \varepsilon_{x} \bX) \bigg]. \end{aligned}

Here, the performance between the recursive SF method and the original SF method can be compared by counting the number of flux function (F(⋅) \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \bF(\cdot)) calls. For instance, the original SF method requires 28 times flux function calls for approximating ∂3F∂U3⋅V⋅W⋅X \newcommand{\bF}{\mathbf{F}} \newcommand{\bU}{\mathbf{U}} \newcommand{\bV}{\mathbf{V}} \newcommand{\bW}{\mathbf{W}} \newcommand{\bX}{\mathbf{X}} \newcommand{\scolon}{\,;\,} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pdd}[2]{\frac{\partial^{2} #1}{\partial #2^{2}}} \newcommand{\pddd}[2]{\frac{\partial^{3} #1}{\partial #2^{3}}} \pddd{\bF}{\bU} \cdot \bV \cdot \bW \cdot \bX term. On the other hand, the recursive system free method needs only eight evaluations. This is a huge improvement in both performance and compactness.

The recursive SF method successfully extended the PIF method [3] to the fourth-order accuracy [1], and you can find the implementations in the SlugCode. The detailed numerical test results are presented in the published article.

[1]Lee, Y., Lee, D., & Reyes, A. (2021). A recursive system-free single-step temporal discretization method for finite difference methods. Journal of Computational Physics: X, 12, 100098. https://doi.org/10.1016/j.jcpx.2021.100098
[2]Lee, Y., & Lee, D. (2021). A single-step third-order temporal discretization with Jacobian-free and Hessian-free formulations for finite difference methods. Journal of Computational Physics, 427, 110063. https://doi.org/10.1016/j.jcp.2020.110063
[3](1, 2) Christlieb, A. J., Guclu, Y., & Seal, D. C. (2015). The Picard integral formulation of weighted essentially nonoscillatory schemes. SIAM Journal on Numerical Analysis, 53(4), 1833-1856. https://doi.org/10.1137/140959936