Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
#pragma once
/**
* Routines for finite differentiation, useful for testing derivatives obtained by other methods
*
* From:
* Bengt Fornberg, 1988, "Coefficients from Generation of Finite Difference Formulas on Arbitrarily Spaced Grids", MATHEMATICS OF COMPUTATION, v. 51, n. 184, pp. 699-706
*
* Higher derivatives should always be done in extended precision mathematics!
*
* Warning: these routines may give entirely erroneous results for double precision arithmetic,
* especially the higher derivatives
*
* Warning: these routines are optimized for accuracy, not for speed or memory use
*/
template<int Nderiv, int Norder, typename Function, typename Scalar>
auto centered_diff(const Function &f, const Scalar x, const Scalar h) {
struct DiffCoeffs {
std::valarray<int> k;
std::valarray<Scalar> c;
};
// See https://en.wikipedia.org/wiki/Finite_difference_coefficient#Central_finite_difference
// Watch out that if you would like to use extended precision, you also need to keep
// all the coefficients in extended precision too
//
// This is because 1.0/2.0 is different than casting each of 1.0 and 2.0 to extended precision
// and then taking their ratio
using r = Scalar;
static std::map<std::tuple<int, int>, DiffCoeffs> CentralDiffCoeffs = {
{{1, 2}, {{-1,1}, {-r(1)/r(2), r(1)/r(2)}} },
{{1, 4}, {{-2,-1,1,2}, {r(1)/ r(1), -r(2)/r(3), r(2)/r(3), -r(1)/r(12)}} },
{{1, 6}, {{-3,-2,-1,1,2,3}, {-r(1)/r(60), r(3)/r(20), -r(3)/r(4), r(3)/r(4), -r(3)/r(20), r(1)/r(60)}} },
{{2, 2}, {{-1,0,1}, {1, -2, 1} }},
{{2, 4}, {{-2,-1,0,1,2}, {r(-1)/r(12), r(4)/r(3), r(-5)/r(2), r(4)/r(3), r(-1)/r(12)}} },
{{2, 6}, {{-3,-2,-1,0,1,2,3}, {r(1)/r(90), r(-3)/r(20), r(3)/r(2), r(-49)/r(18), r(3)/r(2), r(-3)/r(20), r(1)/r(90)}} },
{{3, 2}, {{-2, -1, 0, 1, 2}, {r(-1)/r(2), r(1), 0, r(-1), r(1)/r(2)}} },
{{3, 4}, {{-3,-2,-1,0,1,2,3}, {r(1)/r(8), r(-1), r(13)/r(8), 0, r(-13)/r(8), r(1), r(-1)/r(8)}} },
{{3, 6}, {{-4,-3,-2,-1,0,1,2,3,4}, {r(-7)/r(240), r(3)/r(10), r(-169)/r(120), r(61)/r(30), 0, r(-61)/r(30), r(169)/r(120), r(-3)/r(10), r(7)/r(240)}} },
{{4, 2}, {{-2,-1,0,1,2}, {1,-4,6,-4,1}} },
{{4, 4}, {{-3,-2,-1,0,1,2,3}, {-r(1)/r(6), r(2), -r(13)/r(2), r(28)/r(3.0), -r(13)/r(2), r(2), -r(1)/r(6)}} },
{{4, 6}, {{-4,-3,-2,-1,0,1,2,3,4}, {r(7)/r(240), -r(2)/r(5), r(169)/r(60), -r(122)/r(15), r(91)/r(8), -r(122)/r(15), r(169)/r(60), -r(2)/r(5), r(7)/r(240)}} },
};
auto [k, c] = CentralDiffCoeffs[std::make_tuple(Nderiv, Norder)];
if (c.size() == 0) {
throw std::invalid_argument("Cannot obtain the necessary finite differentiation coefficients");
}
// Sanity check...
if (c.size() != k.size()) {
throw std::invalid_argument("Finite differentiation coefficient arrays not the same size");
}
Scalar num = 0.0;
for (auto i = 0; i < k.size(); ++i) {
num = num + c[i]*f(x + h*k[i]);
}
auto val = num / pow(h, Nderiv);
return val;
}