Add header with finite difference formulae

include/teqp/finite_derivs.hpp

+#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;
+}
\ No newline at end of file