From 9ef68af4677b8e19b9b4dec75ec0c4f086b71f5f Mon Sep 17 00:00:00 2001
From: Timo Koch <timokoch@uio.no>
Date: Wed, 25 Sep 2024 15:49:25 +0200
Subject: [PATCH] [common] Add option for five-point stencil numeric difference
---
dumux/common/numericdifferentiation.hh | 18 +++++++++++++++++-
1 file changed, 17 insertions(+), 1 deletion(-)
diff --git a/dumux/common/numericdifferentiation.hh b/dumux/common/numericdifferentiation.hh
index 9e7af2a78f..7d28a0c46b 100644
--- a/dumux/common/numericdifferentiation.hh
+++ b/dumux/common/numericdifferentiation.hh
@@ -61,7 +61,7 @@ public:
* \param fx0 The result of the function evaluated at x0
* \param eps The numeric epsilon used in the differentiation
* \param numericDifferenceMethod The numeric difference method
- * (1: forward differences (default), 0: central differences, -1: backward differences)
+ * (1: forward differences (default), 0: central differences, -1: backward differences, 5: five-point stencil method)
*/
template<class Function, class Scalar, class FunctionEvalType>
static void partialDerivative(const Function& function, Scalar x0,
@@ -70,7 +70,23 @@ public:
const Scalar eps,
const int numericDifferenceMethod = 1)
{
+ // Five-point stencil numeric difference,
+ // Abramowitz & Stegun, Table 25.2.
+ // The error is proportional to eps^4.
+ if (numericDifferenceMethod == 5)
+ {
+ derivative = function(x0 + eps);
+ derivative -= function(x0 - eps);
+ derivative *= 8.0;
+ derivative += function(x0 - 2*eps);
+ derivative -= function(x0 + 2*eps);
+ derivative /= 12*eps;
+ return;
+ }
+
+ // Forward, central, or backward differences
Scalar delta = 0.0;
+
// we are using forward or central differences, i.e. we need to calculate f(x + \epsilon)
if (numericDifferenceMethod >= 0)
{
--
GitLab