commit 0c4e3e0eb13663e6fcf47211c27096c176701e43
parent 70ebbd4523e821b7b944a37a4de682b81688e8cc
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Sun, 10 Mar 2024 11:06:30 +0100
Updated H calculation
Addition of an internal function H3d_eval that not only calculates H but
also evaluates its first, second and third derivatives. The swf_H3d_eval
function now relies on this function. The derivatives will be used to
evaluate interpolation errors and implement interpolations that are more
efficient than linear interpolation.
Diffstat:
| M | src/swf_H.c | | | 85 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------- |
1 file changed, 70 insertions(+), 15 deletions(-)
diff --git a/src/swf_H.c b/src/swf_H.c
@@ -48,32 +48,87 @@ check_swf_H_tabulate_args(const struct swf_H_tabulate_args* args)
return RES_OK;
}
-/*******************************************************************************
- * Exported symbols
- ******************************************************************************/
- /* H(x) = 1 + 2 Sum(k=1..INF)((-1)^k * exp(-(PI*k)^2 * x)) */
-double
-swf_H3d_eval(const double x)
+/* H(x) = 1 + 2 Sum(k=1..INF)[(-1)^k * exp(-(PI*k)^2 * x)]
+ * H'(x) = - 2 Sum(k=1..INF)[(-1)^k * exp(-(PI*k)^2 * x) * (PI*k)^2]
+ * H''(x) = + 2 Sum(k=1..INF)[(-1)^k * exp(-(PI*k)^2 * x) * (PI*k)^4]
+ * H'''(x) = - 2 Sum(k=1..INF)[(-1)^k * exp(-(PI*k)^2 * x) * (PI*k)^6] */
+static INLINE double
+H3d
+ (const double x,
+ double* dHx, /* Derivative. NULL <=> do not compute it */
+ double* d2Hx, /* Second derivative. NULL <=> do not compute it */
+ double* d3Hx) /* Third derivative. NULL <=> do not compute it */
{
- double sign = -1; /* Sign of the term */
- double sum = 0; /* Sum */
- double error = INF; /* Error */
- unsigned k = 1; /* Index of the term */
+ double Hx = 0; /* Value of the function */
+
+ /* Sum of the terms in the series */
+ double sum = 0; /* Sum of the terms */
+ double sumd = 0; /* Sum of the derivative terms */
+ double sumd2 = 0; /* Sum of the second derivative terms */
+ double sumd3 = 0; /* Sum of the third derivative terms */
+
+ double sign = -1; /* Sign of the serie term */
+ unsigned k = 1; /* Index of the serie term */
+
+ /* Have the calculations converged or are there errors? */
+ int error = 1;
+ int errord = dHx != NULL;
+ int errord2 = d2Hx != NULL;
+ int errord3 = d3Hx != NULL;
+
ASSERT(x >= 0); /* Check pre-condition */
+ /* Do not attempt to calculate 0; it would be numerically catastrophic.
+ * Simply return its value */
if(x == 0) return 0;
do {
const double t = PI*(double)k;
- const double term = sign * exp(-(t*t)*x);
+ const double t2 = t*t;
+ const double term = sign * exp(-t2*x);
const double sum_next = sum + term;
- error = sum_next - sum;
+ error = sum_next != sum;
sum = sum_next;
- sign *= -1;
- } while(error != 0 && ++k < UINT_MAX);
+ sign = -sign;
+
+ if(dHx) { /* Derivative */
+ const double termd = term * t2;
+ const double sumd_next = sumd + termd;
+ errord = sumd_next != sumd;
+ sumd = sumd_next;
+ }
+
+ if(d2Hx) { /* Second derivative */
+ const double termd2 = term * t2*t2;
+ const double sumd2_next = sumd2 + termd2;
+ errord2 = sumd2_next != sumd2;
+ sumd2 = sumd2_next;
+ }
+
+ if(d3Hx) { /* Third derivative */
+ const double termd3 = term * t2*t2*t2;
+ const double sumd3_next = sumd3 + termd3;
+ errord3 = sumd3_next != sumd3;
+ sumd3 = sumd3_next;
+ }
+ } while((error || errord || errord2 || errord3) && ++k < UINT_MAX);
+
+ Hx = 1.0 + 2.0 * sum;
+
+ if(dHx) *dHx = -2.0 * sumd; /* Derivative */
+ if(d2Hx) *d2Hx = 2.0 * sumd2; /* Second derivative */
+ if(d3Hx) *d3Hx = -2.0 * sumd3; /* Third derivative */
+ return Hx;
+}
- return 1.0 + 2.0 * sum;
+/*******************************************************************************
+ * Exported symbols
+ ******************************************************************************/
+double
+swf_H3d_eval(const double x)
+{
+ return H3d(x, NULL, NULL, NULL);
}
double
