fd_coefficients.h
2.25 KB
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
65
66
67
#ifndef STIM_FINITE_DIFFERENCE_COEFFICIENTS
#define STIM_FINITE_DIFFERENCE_COEFFICIENTS
#include <vector>
#include <algorithm>
/// Calculate the coefficients used for finite differences on arbitrary grids. This code is based on:
/// Fornberg, "Generation of Finite Difference Formulas on Arbitrarily Spaced Grids", Mathematics of Computation, 51(184) 699-706, 1988
/// @param Z is the location where the approximation is to be accurate
/// @param X is a vector of coordinates for grid points
/// @param M is the highest derivative to be computed
/// This function returns a matrix C[a][b] where a is the derivative and b is the coefficient for x = b
template <typename T>
std::vector< std::vector<T> > diff_coefficient_matrix(T Z, std::vector<T> X, size_t M){
size_t n = X.size(); //calculate the number of data points available
std::vector< std::vector<T> > C(M+1); //create an array of m arrays
for(size_t m = 0; m < M+1; m++) //for each derivative
C[m].resize(n); //allocate space for n coefficients
T c1, c2, c3, c4, c5;
c1 = 1; //initialize the matrix
c4 = X[0] - Z;
C[0][0] = 1;
for(size_t i = 2; i <= n; i++){
size_t mn = (std::min)(i, M+1);
c2 = 1;
c5 = c4;
c4 = X[i-1] - Z;
for(size_t j = 1; j <= i-1; j++){
c3 = X[i-1] - X[j-1];
c2 = c2 * c3;
if(j == i-1){
for(size_t k = 2; k <= mn; k++){
C[k-1][i-1] = c1 * ((k-1) * C[k-2][i-2] - c5 * C[k-1][i-2]) / c2;
}
C[0][i-1] = -c1 * c5 * C[0][i-2]/c2;
}
T c_p0 = C[0][j-1];
T c_p1 = C[1][j-1];
for(size_t k = 2; k <= mn; k++){
c_p1 = C[k-1][j-1];
C[k-1][j-1] = (c4 * c_p1 - (k-1) * c_p0) / c3;
c_p0 = c_p1;
}
C[0][j-1] = c4 * C[0][j-1] / c3;
}
c1 = c2;
}
return C;
}
/// Returns the coefficients used to estimate the nth derivative at x0 given values at points in x
/// @param x0 is the point where the derivative f'(x0) is to be estimated
/// @param x is the set of points where the values of f(x) are known
/// @param n specifies the derivative, ex. 2 estimates the second derivative f''(x0)
template <typename T>
std::vector<T> diff_coefficients(T x0, std::vector<T> x, size_t n){
std::vector< std::vector<T> > C = diff_coefficient_matrix(x0, x, n); //calculate all derivative coefficients
return C.back();
}
#endif