The system-free (SF) method (see this article) demonstrated its efficiency
and high-order accuracy by combining it with the Picard integration formulation (PIF) method.
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., ),
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,
for the tensor contractions between the same vector, .
In order to approximate the case of different vectors, (, , and for example),
the above equation should be computed multiple times:
As shown above, the original SF method requires
to repeat the approximation seven times
just for getting a single term, .
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 as,
which is the Jacobian-free approximation from the original SF method.
For higher-order derivatives, we take the functional
in a successive fashion. For example, the Hessian-vector contraction would be,
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 .
Following the equivalent strategy,
the tensor contraction of the third-order derivative of the flux function,
with three distinct vectors, and
can be expressed in a compact form as,
Here, the performance between the recursive SF method and the original SF method
can be compared by counting the number of flux function () calls.
For instance, the original SF method requires 28 times flux function calls
for approximating 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 to the fourth-order accuracy ,
and you can find the implementations in the SlugCode.
The detailed numerical test results are presented in the published article.