#ifndef STIM_FINITE_DIFFERENCE_COEFFICIENTS #define STIM_FINITE_DIFFERENCE_COEFFICIENTS #include #include /// 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 std::vector< std::vector > diff_coefficient_matrix(T Z, std::vector X, size_t M){ size_t n = X.size(); //calculate the number of data points available std::vector< std::vector > 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 std::vector diff_coefficients(T x0, std::vector x, size_t n){ std::vector< std::vector > C = diff_coefficient_matrix(x0, x, n); //calculate all derivative coefficients return C.back(); } #endif