How can I estimate the first derivative of an evaluated univariate function in C++?

3
votes

I have a two vectors:

  • Coordinates along an axis, x;
  • An evaluation of a function on those coordinates, f(x).

and I want to compute an estimate of the first derivative of f at these coordinates.

The function is a descriptor of a wavefunction and x is the dihedral angle.

Because the result vector must have the same length as the two existing vectors, I cannot use a manual implementation based on Newton's difference quotient.

In Python, I can obtain such an estimate using the scipy library:

spline = UnivariateSpline(X, Y, k=4, s=0)
sd = spline.derivative(n=1)

Perhaps I can do something similar in C++?

1
c++numerical-methodsderivative
Please don't describe your code, create a minimal reproducible example of it to show us instead. And please take some time to read or refresh the help pages, take the SO tour, read How to Ask, as well as this question checklist. - Some programmer dude
There is nothing for this in the C++ standard library. Mathematically speaking, what is x and what is f(x)? Do you know what the derivative of f with respect to x is if x and f(x) are scalars? If x is a scalar and f(x) is a vector? If x is a vector and f(x) is a scalar? If both are vectors? Know what you want conceptually before throwing code at it. - gspr
Both f(x) and x are scalars and I don't know the exact form of f(x), but only its values. - Jonny_92
@TheOldJonny: Your Python code uses a degree-4 spline, which it creates based on the knots. You should search - using a search engine, and perhaps also on GitHub - for C/C++ libraries for working with such splines. The Eigen library has a spline module which may fit the bill (not sure - haven't used it). - einpoklum
If you want to write your own code, you may use the formulae for numerical differentiation, since you have also the xvector, but you have to consider some "boundary conditions". The solution you wrote in python consists of creating a fit of the original function and then deriving the latter, you are not computing the derivatives of the original function (although you may achieve a good rate of approximation depending on f). - Eddymage

1 Answers

2
votes

I'd start with pchip.

Example:

#include <boost/math/interpolators/pchip.hpp>
// ...
using boost::math::interpolators::pchip;
auto f = pchip(std::move(x), std::move(y));
double t = 3.2;
std::cout << "f(" << t << " = " << f(t) << ", f'(" << t << ") = " << f.prime(t) << "\n";

If you don't like the "character" of pchip, then you have many other options.