fd_coefficients.h 2.25 KB
#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