c ++ Eigen：样条导数（）给出奇怪的导数
c++ Eigen: spline derivatives() gives strange derivatives
我正在尝试了解如何在 Eigen 中使用样条曲线，特别是我想在某个点上找到样条插值及其一阶和二阶导数的值。找到内插值很容易，但是当我尝试计算导数时，我得到了奇怪的值。
我尝试按照手册中 derivatives 命令的说明进行操作 (http://eigen.tuxfamily.org/dox/unsupported/classEigen_1_1Spline.html#af3586ab1929959e0161bfe7da40155c6)，这是我在代码中的尝试：
#include <iostream>
#include <Eigen/Core>
#include <unsupported/Eigen/Splines>
using namespace Eigen;
using namespace std;
double scaling(double x, double min, double max) // for scaling numbers
{
return (x - min)/(max - min);
}
VectorXd scale(VectorXd xvals) // for scaling vectors
{
const double min = xvals.minCoeff();
const double max = xvals.maxCoeff();
for (int k = 0; k < xvals.size(); k++)
xvals(k) = scaling(xvals(k),min,max);
return xvals;
}
int main()
{
typedef Spline<double,1,3> spline;
VectorXd xvals = (VectorXd(4) << 0,1,2,4).finished();
VectorXd yvals = xvals.array().square(); // x^2
spline testspline = SplineFitting<spline>::Interpolate(yvals.transpose(), 3,
scale(xvals).transpose());
cout << "derivative at x = 0: " << testspline.derivatives(0.00,2) << endl;
cout << "derivative at x = 1: " << testspline.derivatives(0.25,2) << endl;
cout << "derivative at x = 2: " << testspline.derivatives(0.50,2) << endl;
cout << "derivative at x = 3: " << testspline.derivatives(0.75,2) << endl;
cout << "derivative at x = 4: " << testspline.derivatives(1.00,2) << endl;
}
它输出
derivative at x = 0: 0 0 32
derivative at x = 1: 1 8 32
derivative at x = 2: 4 16 32
derivative at x = 3: 9 24 32
derivative at x = 4: 16 32 32
也就是插值是正确的(c.f.x = 3)，但是导数不正确，而且是系统性的偏离，所以我在想我做错了什么.由于这些遵循 x^2，导数应为 0,2,4,6,8，二阶导数应为 2.
关于如何解决这个问题有什么想法吗？
编辑 1
将 x^2 更改为 x^2 + 1 会产生相同的导数，因此至少可以检查出来。但是将 x^2 更改为 x^3 是错误的，但错误的方式略有不同，输出将是：
derivative at x = 2: 8 48 192
derivative at x = 3: 27 108 288
derivative at x = 4: 64 192 384
错了，应该是6, 9, 12。
还有运行 x^2 的情况，但是将输入向量更改为 0,1,...9 会产生与使用原始输入向量相同的导数，但二阶导数变为稳定的 200，这也是错误的。我不明白为什么二阶导数应该取决于输入点的数量。
编辑：参见工作 .h-file 计算底部任意订单的 B-splines。
免责声明：这不是我的问题的答案，因为它实际上已在标题中说明，而是 work-around 和一些评论。
经过与用户@Paul H. 的商议（见评论），我意识到我对样条曲线的有限理解可能对我造成了一些困惑。在对 Eigen 文档进行一些审查之后， derivative() 命令确实按预期工作似乎是合理的，因此使我的问题措辞不当。 derivative() 计算 样条曲线 的导数，而不是像我预期的那样计算拟合 函数 的导数。我还没有想出一种方法让 Eigen 从拟合中输出函数导数，我认为它不是为此而设计的。然而，一旦使用一些计算导数的标准算法获得拟合点，当然可以很容易地计算导数。
我写了下面的.h-file用于计算样条及其导数的过程，我认为值得分享。为了方便起见，它的评论相当好。
请注意，此程序使用样条的 1 索引而不是 0 索引，因此对于例如二次 B-spline order 应设置为 4。这是一个小问题，可以通过更改计算以匹配 wikipedia.
轻松修复
bsplines.h
//
// Header-file for calculating splines using standard libraries
//
// usage:
//
// x bsplines class constructs a set of splines up to some given order,
// on some given knot sequence, the splines are stored in a vector,
// such that splines[a][b] accesses the spline of order a and index b
// x get<some_member>() is an accessor that returns a pointer to some
// data member of the spline
// x calcsplines() calculates spline values as well as first and second
// order derivatives on some predefined grid
// x calcspline() returns the spline value as well as first and second
// derivatives in some point. This alborithm is slower than the grid
// one, due to unnecessary recalculations of intermediate results
// x writesplines() writes the splines and their derivatives to a file
// x for more details se the class declaration below
// TODO:
// x change to 0-indexation
// x introduce the possibility of calculating higher order derivatives
// recursively
//
// change log:
//
// 1.0 - initial release
// 1.1 - reworked grid such that the class now expects separate
// grid and knot files.
// - added the ability to calculate spline value in a point
// rather than calculate values on a grid
// - added a feature to change knots and grid
// 1.1.1 - reworked how returning single values works
// 1.1.2 - enabled swapping grid
//
// Note:
//
// This file uses 1-indexation rathar than 0-indexation, hence a qubic spline
// would be k = 4. Someone should eventually fix this as this is non-standard.
//
// Also, while only standard libraries are used here, you might want to check out
// some linear algebra package (e.g. Armadillo or Eigen) if you're going to use the
// splines in a context where you need linear algebraic operations.
//
// Originally developed by David Andersson
//
#include <iomanip>
#include <sstream>
#include <iostream>
#include <vector>
#include <algorithm>
#include <fstream>
#include <functional>
using namespace std;
typedef unsigned int uint;
class bsplines // class for bsplines
{
// private section
uint order; // order of spline
uint gridpts; // number of grid points
uint knotpts; // number of knot points
double tolerance; // tolerance for float comparisons
vector<double> knots; // knot sequence
vector<double> grid; // grid points
class spline // a member spline in the set of splines
{
int index; // the spline index, or number
vector<double> vals; // spline values
vector<double> d1; // spline first derivatives
vector<double> d2; // spline second derivatives
double tval; // same, but in one point
double td1;
double td2;
friend bsplines; // for ease of access
public:
};
vector<vector <spline>> splines; // the set of splines
// puclic section
public:
void readknots(string); // read knots from file
void readknotsnorm(string); // read knots from file and normalize
void readgrid(string); // read grid from file
void swapgrid(string); // reads and swaps new grid from file
void writesplines(); // write spline vals and derivs to file
void buildsplines(); // build the set of splines
void calcsplines(); // calculate spline vals and derivs
void printknots(); // print knot sequence
void printgrid(); // print grid
void printgridsize(); // print gridsize
void printvals(uint,uint); // print values of a spline
vector <double> calcspline(uint,uint,double); // calculate spline in point
// accessors // returns pointer to member
vector <double>* getknots(){return &knots;}
vector <double>* getgrid(){return &grid;}
uint* getknotpts(){return &knotpts;}
uint* getgridpts(){return &gridpts;}
uint getnosplines(uint m){return splines[m].size();}
vector <spline>* getsplines(uint m){return &splines[m];}
vector <double>* getvals(uint m, uint n){return &splines[m][n].vals;}
vector <double>* getd1(uint m, uint n){return &splines[m][n].d1;}
vector <double>* getd2(uint m, uint n){return &splines[m][n].d2;}
// constructor // sets up the spline class
bsplines (string iknots, string igrid, uint iorder, double itol)
:order(iorder), tolerance(itol)
{
readknots(iknots);
readgrid(igrid);
buildsplines();
}
};
void bsplines::buildsplines()
{
{
for (uint l = 1; l <= order; l++)
{
vector <spline> splinevec;
for (uint k = 0; k < knotpts - l; k++)
{
spline tmp;
tmp.index = k;
tmp.vals.reserve(gridpts);
tmp.d1.reserve(gridpts);
tmp.d2.reserve(gridpts);
splinevec.push_back(tmp);
}
splines.push_back(splinevec);
}
}
}
vector <double> bsplines::calcspline(uint m, uint n, double x)
{
// first order splines // exceptions handles infinities
for (auto& sp : splines[0])
{
uint i = sp.index;
if (x > knots[i+1])
sp.tval = 0;
else if ((x >= knots[i] && x < knots[i+1]) || x == knots.back())
sp.tval = 1;
else
sp.tval = 0;
}
// higher order splines
for (uint o = 1; o < order; o++)
{
uint oo = o+1; // compensating for 1-indexation
for (auto& sp : splines[o])
{
uint i = sp.index;
double t1 = knots[i+oo-1] - knots[i];
double t2 = knots[i+oo] - knots[i+1];
double c = 0;
if (abs(t1) > tolerance)
c += (x - knots[i]) / t1 * splines[o-1][i].tval;
if (abs(t2) > tolerance)
c += (knots[i+oo] - x) / t2 * splines[o-1][i+1].tval;
sp.tval = c;
}
}
uint o = order - 1;
// first order derivatives
for (auto& sp : splines[o])
{
uint i = sp.index;
double t1 = knots[i+order-1] - knots[i];
double t2 = knots[i+order] - knots[i+1];
double c = 0;
if (abs(t1) > tolerance)
c += 1.0 / t1 * splines[o-1][i].tval;
if (abs(t2) > tolerance)
c -= 1.0 / t2 * splines[o-1][i+1].tval;
c *= (order-1);
sp.td1 = c;
}
// second order derivatives
for (auto& sp : splines[o])
{
uint i = sp.index;
double t1 = (knots[i+order-1] - knots[i+0]) * (knots[i+order-2] - knots[i+0]);
double t2 = (knots[i+order-1] - knots[i+0]) * (knots[i+order-1] - knots[i+1]);
double t3 = (knots[i+order-0] - knots[i+1]) * (knots[i+order-1] - knots[i+1]);
double t4 = (knots[i+order-0] - knots[i+1]) * (knots[i+order-0] - knots[i+2]);
double c = 0;
if (abs(t1) > tolerance)
c += 1.0 / t1 * splines[o-2][sp.index].tval;
if (abs(t2) > tolerance)
c -= 1.0 / t2 * splines[o-2][sp.index+1].tval;
if (abs(t3) > tolerance)
c -= 1.0 / t3 * splines[o-2][sp.index+1].tval;
if (abs(t4) > tolerance)
c += 1.0 / t4 * splines[o-2][sp.index+2].tval;
c *= (order-1)*(order-2);
sp.td2 = c;
}
vector <double> retvals = {splines[m][n].tval, splines[m][n].td1, splines[m][n].td2};
return retvals;
}
void bsplines::calcsplines()
{
// first order splines
for (auto& sp : splines[0])
{
uint i = sp.index;
for (auto& x : grid)
{
if (x > knots[i+1])
sp.vals.push_back(0);
else if ((x >= knots[i] && x < knots[i+1]) || x == knots.back())
sp.vals.push_back(1);
else
sp.vals.push_back(0);
}
}
// higher order splines
for (uint o = 1; o < order; o++)
{
uint oo = o+1; // compensating for 1-indexation
for (auto& sp : splines[o])
{
uint i = sp.index;
double t1 = knots[i+oo-1] - knots[i];
double t2 = knots[i+oo] - knots[i+1];
for (auto& x : grid)
{
uint k = &x - &grid[0];
double c = 0;
if (abs(t1) > tolerance)
c += (x - knots[i]) / t1 * splines[o-1][i].vals[k];
if (abs(t2) > tolerance)
c += (knots[i+oo] - x) / t2 * splines[o-1][i+1].vals[k];
sp.vals.push_back(c);
}
}
}
uint o = order - 1; // use this one when accessing splines;
// first order derivatives
for (auto& sp : splines[o])
{
uint i = sp.index;
double t1 = knots[i+order-1] - knots[i];
double t2 = knots[i+order] - knots[i+1];
for (auto& x : grid)
{
uint k = &x - &grid[0];
double c = 0;
if (abs(t1) > tolerance)
c += 1.0 / t1 * splines[o-1][i].vals[k];
if (abs(t2) > tolerance)
c -= 1.0 / t2 * splines[o-1][i+1].vals[k];
c *= (order-1);
sp.d1.push_back(c);
}
}
// second order derivatives
for (auto& sp : splines[o])
{
uint i = sp.index;
double t1 = (knots[i+order-1] - knots[i+0]) * (knots[i+order-2] - knots[i+0]);
double t2 = (knots[i+order-1] - knots[i+0]) * (knots[i+order-1] - knots[i+1]);
double t3 = (knots[i+order-0] - knots[i+1]) * (knots[i+order-1] - knots[i+1]);
double t4 = (knots[i+order-0] - knots[i+1]) * (knots[i+order-0] - knots[i+2]);
for (auto& x : grid)
{
uint k = &x - &grid[0];
double c = 0;
if (abs(t1) > tolerance)
c += 1.0 / t1 * splines[o-2][sp.index].vals[k];
if (abs(t2) > tolerance)
c -= 1.0 / t2 * splines[o-2][sp.index+1].vals[k];
if (abs(t3) > tolerance)
c -= 1.0 / t3 * splines[o-2][sp.index+1].vals[k];
if (abs(t4) > tolerance)
c += 1.0 / t4 * splines[o-2][sp.index+2].vals[k];
c *= (order-1)*(order-2);
sp.d2.push_back(c);
}
}
}
void bsplines::readknots(string knotfile)
{
double x;
ifstream readknots(knotfile);
while (readknots >> x)
knots.push_back(x);
for (uint k = 0; k < order - 1; k++)
{
knots.insert(knots.begin(),knots.front());
knots.insert(knots.end(),knots.back());
}
knotpts = knots.size();
}
void bsplines::readknotsnorm(string knotfile)
{
double x;
knots.reserve(knotpts + 2*(order - 1));
ifstream readknots(knotfile);
while (readknots >> x)
knots.push_back(x);
auto minmax = minmax_element(begin(knots), end(knots));
double min = *(minmax.first);
double max = *(minmax.second);
for (auto& el : knots)
el = (el - min) / (max-min);
}
void bsplines::readgrid(string gridfile)
{
double x;
ifstream readgrid(gridfile);
while (readgrid >> x)
grid.push_back(x);
gridpts = grid.size();
}
void bsplines::swapgrid(string gridfile)
{
grid = {};
double x;
ifstream readgrid(gridfile);
while (readgrid >> x)
grid.push_back(x);
gridpts = grid.size();
}
void bsplines::printknots()
{
cout << "content in knot vector: " << endl;
for (auto& el : knots)
cout << el << " ";
cout << endl;
}
void bsplines::printgrid()
{
cout << "content in grid vector: " << endl;
for (auto& el : grid)
cout << el << " ";
cout << endl;
}
void bsplines::printgridsize()
{
cout << "number of grid points: " << endl << grid.size() << endl;
}
void bsplines::printvals(uint m, uint n)
{
cout << "content in spline (B" << m << "," << n << ") vals vector: " << endl;
for (auto& el : splines[n][m].vals)
cout << el << " ";
cout << endl;
}
void bsplines::writesplines()
{
for (uint o = 0; o < order; o++)
for (auto& sp : splines[o])
{
uint i = sp.index;
ostringstream namestream;
namestream << "B(" << fixed << setprecision(1) << i << ","
<< fixed << setprecision(1) << o << ").csv";
string filename = namestream.str();
ofstream fs;
fs.open(filename);
if (o < order - 1)
{
for (uint k = 0; k < sp.vals.size(); k++)
fs << sp.vals[k] << "," << 0 << "," << 0 << endl;
fs.close();
}
else
{
for (uint k = 0; k < sp.vals.size(); k++)
fs << sp.vals[k] << "," << sp.d1[k] << "," << sp.d2[k] << endl;
fs.close();
}
cout << "write " << sp.vals.size() << " numbers to " << filename << endl;
}
}
编辑：已更新。h-file。
解决了。你们非常亲密。您所要做的就是扩展衍生品
与
1 / (x_max - x_min)（一阶导数）1 / (x_max - x_min)^2（二阶导数）。
TLDR：您在拟合样条曲线时将 x 值归一化为介于 0 和 1 之间，但没有缩放 y 值。
您实际拟合的不是样条拟合 x^2，而是：
x_norm = (x - x_min) / (x_max - x_min)
y = x_norm**2
因此使用链式法则 y = x_norm**2 的一阶导数为 2x / (x_max - x_min)，二阶导数为 2 / (x_max - x_min)**2.
完整示例代码：
#include <iostream>
#include <Eigen/Core>
#include <unsupported/Eigen/Splines>
using namespace Eigen;
using namespace std;
VectorXd normalize(const VectorXd &x) {
VectorXd x_norm;
x_norm.resize(x.size());
const double min = x.minCoeff();
const double max = x.maxCoeff();
for (int k = 0; k < x.size(); k++) {
x_norm(k) = (x(k) - min)/(max - min);
}
return x_norm;
}
int main() {
typedef Spline<double, 1, 3> Spline1D;
typedef SplineFitting<Spline1D> Spline1DFitting;
const Vector4d x{0, 1, 2, 4};
const Vector4d y = (x.array().square()); // x^2
const auto knots = normalize(x); // Normalize x to be between 0 and 1
const double scale = 1 / (x.maxCoeff() - x.minCoeff());
const double scale_sq = scale * scale;
Spline1D spline = Spline1DFitting::Interpolate(y.transpose(), 3, knots);
cout << "1st deriv at x = 0: " << spline.derivatives(0.00, 1)(1) * scale << endl;
cout << "1st deriv at x = 1: " << spline.derivatives(0.25, 1)(1) * scale << endl;
cout << "1st deriv at x = 2: " << spline.derivatives(0.50, 1)(1) * scale << endl;
cout << "1st deriv at x = 3: " << spline.derivatives(0.75, 1)(1) * scale << endl;
cout << "1st deriv at x = 4: " << spline.derivatives(1.00, 1)(1) * scale << endl;
cout << endl;
/**
* IMPORTANT NOTE: Eigen's spline module is not documented well. Once you fit a spline
* to find the derivative of the fitted spline at any point u [0, 1] you call:
*
* spline.derivatives(u, 1)(1)
* ^ ^ ^
* | | |
* | | +------- Access the result
* | +---------- Derivative order
* +------------- Parameter u [0, 1]
*
* The last bit `(1)` is if the spline is 1D. And value of `1` for the first
* order. `2` for the second order. Do not forget to scale the result.
*
* For higher dimensions, treat the return as a matrix and grab the 1st or
* 2nd column for the first and second derivative.
*/
cout << "2nd deriv at x = 0: " << spline.derivatives(0.00, 2)(2) * scale_sq << endl;
cout << "2nd deriv at x = 1: " << spline.derivatives(0.25, 2)(2) * scale_sq << endl;
cout << "2nd deriv at x = 2: " << spline.derivatives(0.50, 2)(2) * scale_sq << endl;
cout << "2nd deriv at x = 3: " << spline.derivatives(0.75, 2)(2) * scale_sq << endl;
cout << "2nd deriv at x = 4: " << spline.derivatives(1.00, 2)(2) * scale_sq << endl;
return 0;
}
示例输出：
1st deriv at x = 0: 4.52754e-16
1st deriv at x = 1: 2
1st deriv at x = 2: 4
1st deriv at x = 3: 6
1st deriv at x = 4: 8
2nd deriv at x = 0: 2
2nd deriv at x = 1: 2
2nd deriv at x = 2: 2
2nd deriv at x = 3: 2
2nd deriv at x = 4: 2
我正在尝试了解如何在 Eigen 中使用样条曲线，特别是我想在某个点上找到样条插值及其一阶和二阶导数的值。找到内插值很容易，但是当我尝试计算导数时，我得到了奇怪的值。
我尝试按照手册中 derivatives 命令的说明进行操作 (http://eigen.tuxfamily.org/dox/unsupported/classEigen_1_1Spline.html#af3586ab1929959e0161bfe7da40155c6)，这是我在代码中的尝试：
#include <iostream>
#include <Eigen/Core>
#include <unsupported/Eigen/Splines>
using namespace Eigen;
using namespace std;
double scaling(double x, double min, double max) // for scaling numbers
{
return (x - min)/(max - min);
}
VectorXd scale(VectorXd xvals) // for scaling vectors
{
const double min = xvals.minCoeff();
const double max = xvals.maxCoeff();
for (int k = 0; k < xvals.size(); k++)
xvals(k) = scaling(xvals(k),min,max);
return xvals;
}
int main()
{
typedef Spline<double,1,3> spline;
VectorXd xvals = (VectorXd(4) << 0,1,2,4).finished();
VectorXd yvals = xvals.array().square(); // x^2
spline testspline = SplineFitting<spline>::Interpolate(yvals.transpose(), 3,
scale(xvals).transpose());
cout << "derivative at x = 0: " << testspline.derivatives(0.00,2) << endl;
cout << "derivative at x = 1: " << testspline.derivatives(0.25,2) << endl;
cout << "derivative at x = 2: " << testspline.derivatives(0.50,2) << endl;
cout << "derivative at x = 3: " << testspline.derivatives(0.75,2) << endl;
cout << "derivative at x = 4: " << testspline.derivatives(1.00,2) << endl;
}
它输出
derivative at x = 0: 0 0 32
derivative at x = 1: 1 8 32
derivative at x = 2: 4 16 32
derivative at x = 3: 9 24 32
derivative at x = 4: 16 32 32
也就是插值是正确的(c.f.x = 3)，但是导数不正确，而且是系统性的偏离，所以我在想我做错了什么.由于这些遵循 x^2，导数应为 0,2,4,6,8，二阶导数应为 2.
关于如何解决这个问题有什么想法吗？
编辑 1
将 x^2 更改为 x^2 + 1 会产生相同的导数，因此至少可以检查出来。但是将 x^2 更改为 x^3 是错误的，但错误的方式略有不同，输出将是：
derivative at x = 2: 8 48 192
derivative at x = 3: 27 108 288
derivative at x = 4: 64 192 384
错了，应该是6, 9, 12。
还有运行 x^2 的情况，但是将输入向量更改为 0,1,...9 会产生与使用原始输入向量相同的导数，但二阶导数变为稳定的 200，这也是错误的。我不明白为什么二阶导数应该取决于输入点的数量。
编辑：参见工作 .h-file 计算底部任意订单的 B-splines。
免责声明：这不是我的问题的答案，因为它实际上已在标题中说明，而是 work-around 和一些评论。
经过与用户@Paul H. 的商议（见评论），我意识到我对样条曲线的有限理解可能对我造成了一些困惑。在对 Eigen 文档进行一些审查之后， derivative() 命令确实按预期工作似乎是合理的，因此使我的问题措辞不当。 derivative() 计算 样条曲线 的导数，而不是像我预期的那样计算拟合 函数 的导数。我还没有想出一种方法让 Eigen 从拟合中输出函数导数，我认为它不是为此而设计的。然而，一旦使用一些计算导数的标准算法获得拟合点，当然可以很容易地计算导数。
我写了下面的.h-file用于计算样条及其导数的过程，我认为值得分享。为了方便起见，它的评论相当好。
请注意，此程序使用样条的 1 索引而不是 0 索引，因此对于例如二次 B-spline order 应设置为 4。这是一个小问题，可以通过更改计算以匹配 wikipedia.
bsplines.h
//
// Header-file for calculating splines using standard libraries
//
// usage:
//
// x bsplines class constructs a set of splines up to some given order,
// on some given knot sequence, the splines are stored in a vector,
// such that splines[a][b] accesses the spline of order a and index b
// x get<some_member>() is an accessor that returns a pointer to some
// data member of the spline
// x calcsplines() calculates spline values as well as first and second
// order derivatives on some predefined grid
// x calcspline() returns the spline value as well as first and second
// derivatives in some point. This alborithm is slower than the grid
// one, due to unnecessary recalculations of intermediate results
// x writesplines() writes the splines and their derivatives to a file
// x for more details se the class declaration below
// TODO:
// x change to 0-indexation
// x introduce the possibility of calculating higher order derivatives
// recursively
//
// change log:
//
// 1.0 - initial release
// 1.1 - reworked grid such that the class now expects separate
// grid and knot files.
// - added the ability to calculate spline value in a point
// rather than calculate values on a grid
// - added a feature to change knots and grid
// 1.1.1 - reworked how returning single values works
// 1.1.2 - enabled swapping grid
//
// Note:
//
// This file uses 1-indexation rathar than 0-indexation, hence a qubic spline
// would be k = 4. Someone should eventually fix this as this is non-standard.
//
// Also, while only standard libraries are used here, you might want to check out
// some linear algebra package (e.g. Armadillo or Eigen) if you're going to use the
// splines in a context where you need linear algebraic operations.
//
// Originally developed by David Andersson
//
#include <iomanip>
#include <sstream>
#include <iostream>
#include <vector>
#include <algorithm>
#include <fstream>
#include <functional>
using namespace std;
typedef unsigned int uint;
class bsplines // class for bsplines
{
// private section
uint order; // order of spline
uint gridpts; // number of grid points
uint knotpts; // number of knot points
double tolerance; // tolerance for float comparisons
vector<double> knots; // knot sequence
vector<double> grid; // grid points
class spline // a member spline in the set of splines
{
int index; // the spline index, or number
vector<double> vals; // spline values
vector<double> d1; // spline first derivatives
vector<double> d2; // spline second derivatives
double tval; // same, but in one point
double td1;
double td2;
friend bsplines; // for ease of access
public:
};
vector<vector <spline>> splines; // the set of splines
// puclic section
public:
void readknots(string); // read knots from file
void readknotsnorm(string); // read knots from file and normalize
void readgrid(string); // read grid from file
void swapgrid(string); // reads and swaps new grid from file
void writesplines(); // write spline vals and derivs to file
void buildsplines(); // build the set of splines
void calcsplines(); // calculate spline vals and derivs
void printknots(); // print knot sequence
void printgrid(); // print grid
void printgridsize(); // print gridsize
void printvals(uint,uint); // print values of a spline
vector <double> calcspline(uint,uint,double); // calculate spline in point
// accessors // returns pointer to member
vector <double>* getknots(){return &knots;}
vector <double>* getgrid(){return &grid;}
uint* getknotpts(){return &knotpts;}
uint* getgridpts(){return &gridpts;}
uint getnosplines(uint m){return splines[m].size();}
vector <spline>* getsplines(uint m){return &splines[m];}
vector <double>* getvals(uint m, uint n){return &splines[m][n].vals;}
vector <double>* getd1(uint m, uint n){return &splines[m][n].d1;}
vector <double>* getd2(uint m, uint n){return &splines[m][n].d2;}
// constructor // sets up the spline class
bsplines (string iknots, string igrid, uint iorder, double itol)
:order(iorder), tolerance(itol)
{
readknots(iknots);
readgrid(igrid);
buildsplines();
}
};
void bsplines::buildsplines()
{
{
for (uint l = 1; l <= order; l++)
{
vector <spline> splinevec;
for (uint k = 0; k < knotpts - l; k++)
{
spline tmp;
tmp.index = k;
tmp.vals.reserve(gridpts);
tmp.d1.reserve(gridpts);
tmp.d2.reserve(gridpts);
splinevec.push_back(tmp);
}
splines.push_back(splinevec);
}
}
}
vector <double> bsplines::calcspline(uint m, uint n, double x)
{
// first order splines // exceptions handles infinities
for (auto& sp : splines[0])
{
uint i = sp.index;
if (x > knots[i+1])
sp.tval = 0;
else if ((x >= knots[i] && x < knots[i+1]) || x == knots.back())
sp.tval = 1;
else
sp.tval = 0;
}
// higher order splines
for (uint o = 1; o < order; o++)
{
uint oo = o+1; // compensating for 1-indexation
for (auto& sp : splines[o])
{
uint i = sp.index;
double t1 = knots[i+oo-1] - knots[i];
double t2 = knots[i+oo] - knots[i+1];
double c = 0;
if (abs(t1) > tolerance)
c += (x - knots[i]) / t1 * splines[o-1][i].tval;
if (abs(t2) > tolerance)
c += (knots[i+oo] - x) / t2 * splines[o-1][i+1].tval;
sp.tval = c;
}
}
uint o = order - 1;
// first order derivatives
for (auto& sp : splines[o])
{
uint i = sp.index;
double t1 = knots[i+order-1] - knots[i];
double t2 = knots[i+order] - knots[i+1];
double c = 0;
if (abs(t1) > tolerance)
c += 1.0 / t1 * splines[o-1][i].tval;
if (abs(t2) > tolerance)
c -= 1.0 / t2 * splines[o-1][i+1].tval;
c *= (order-1);
sp.td1 = c;
}
// second order derivatives
for (auto& sp : splines[o])
{
uint i = sp.index;
double t1 = (knots[i+order-1] - knots[i+0]) * (knots[i+order-2] - knots[i+0]);
double t2 = (knots[i+order-1] - knots[i+0]) * (knots[i+order-1] - knots[i+1]);
double t3 = (knots[i+order-0] - knots[i+1]) * (knots[i+order-1] - knots[i+1]);
double t4 = (knots[i+order-0] - knots[i+1]) * (knots[i+order-0] - knots[i+2]);
double c = 0;
if (abs(t1) > tolerance)
c += 1.0 / t1 * splines[o-2][sp.index].tval;
if (abs(t2) > tolerance)
c -= 1.0 / t2 * splines[o-2][sp.index+1].tval;
if (abs(t3) > tolerance)
c -= 1.0 / t3 * splines[o-2][sp.index+1].tval;
if (abs(t4) > tolerance)
c += 1.0 / t4 * splines[o-2][sp.index+2].tval;
c *= (order-1)*(order-2);
sp.td2 = c;
}
vector <double> retvals = {splines[m][n].tval, splines[m][n].td1, splines[m][n].td2};
return retvals;
}
void bsplines::calcsplines()
{
// first order splines
for (auto& sp : splines[0])
{
uint i = sp.index;
for (auto& x : grid)
{
if (x > knots[i+1])
sp.vals.push_back(0);
else if ((x >= knots[i] && x < knots[i+1]) || x == knots.back())
sp.vals.push_back(1);
else
sp.vals.push_back(0);
}
}
// higher order splines
for (uint o = 1; o < order; o++)
{
uint oo = o+1; // compensating for 1-indexation
for (auto& sp : splines[o])
{
uint i = sp.index;
double t1 = knots[i+oo-1] - knots[i];
double t2 = knots[i+oo] - knots[i+1];
for (auto& x : grid)
{
uint k = &x - &grid[0];
double c = 0;
if (abs(t1) > tolerance)
c += (x - knots[i]) / t1 * splines[o-1][i].vals[k];
if (abs(t2) > tolerance)
c += (knots[i+oo] - x) / t2 * splines[o-1][i+1].vals[k];
sp.vals.push_back(c);
}
}
}
uint o = order - 1; // use this one when accessing splines;
// first order derivatives
for (auto& sp : splines[o])
{
uint i = sp.index;
double t1 = knots[i+order-1] - knots[i];
double t2 = knots[i+order] - knots[i+1];
for (auto& x : grid)
{
uint k = &x - &grid[0];
double c = 0;
if (abs(t1) > tolerance)
c += 1.0 / t1 * splines[o-1][i].vals[k];
if (abs(t2) > tolerance)
c -= 1.0 / t2 * splines[o-1][i+1].vals[k];
c *= (order-1);
sp.d1.push_back(c);
}
}
// second order derivatives
for (auto& sp : splines[o])
{
uint i = sp.index;
double t1 = (knots[i+order-1] - knots[i+0]) * (knots[i+order-2] - knots[i+0]);
double t2 = (knots[i+order-1] - knots[i+0]) * (knots[i+order-1] - knots[i+1]);
double t3 = (knots[i+order-0] - knots[i+1]) * (knots[i+order-1] - knots[i+1]);
double t4 = (knots[i+order-0] - knots[i+1]) * (knots[i+order-0] - knots[i+2]);
for (auto& x : grid)
{
uint k = &x - &grid[0];
double c = 0;
if (abs(t1) > tolerance)
c += 1.0 / t1 * splines[o-2][sp.index].vals[k];
if (abs(t2) > tolerance)
c -= 1.0 / t2 * splines[o-2][sp.index+1].vals[k];
if (abs(t3) > tolerance)
c -= 1.0 / t3 * splines[o-2][sp.index+1].vals[k];
if (abs(t4) > tolerance)
c += 1.0 / t4 * splines[o-2][sp.index+2].vals[k];
c *= (order-1)*(order-2);
sp.d2.push_back(c);
}
}
}
void bsplines::readknots(string knotfile)
{
double x;
ifstream readknots(knotfile);
while (readknots >> x)
knots.push_back(x);
for (uint k = 0; k < order - 1; k++)
{
knots.insert(knots.begin(),knots.front());
knots.insert(knots.end(),knots.back());
}
knotpts = knots.size();
}
void bsplines::readknotsnorm(string knotfile)
{
double x;
knots.reserve(knotpts + 2*(order - 1));
ifstream readknots(knotfile);
while (readknots >> x)
knots.push_back(x);
auto minmax = minmax_element(begin(knots), end(knots));
double min = *(minmax.first);
double max = *(minmax.second);
for (auto& el : knots)
el = (el - min) / (max-min);
}
void bsplines::readgrid(string gridfile)
{
double x;
ifstream readgrid(gridfile);
while (readgrid >> x)
grid.push_back(x);
gridpts = grid.size();
}
void bsplines::swapgrid(string gridfile)
{
grid = {};
double x;
ifstream readgrid(gridfile);
while (readgrid >> x)
grid.push_back(x);
gridpts = grid.size();
}
void bsplines::printknots()
{
cout << "content in knot vector: " << endl;
for (auto& el : knots)
cout << el << " ";
cout << endl;
}
void bsplines::printgrid()
{
cout << "content in grid vector: " << endl;
for (auto& el : grid)
cout << el << " ";
cout << endl;
}
void bsplines::printgridsize()
{
cout << "number of grid points: " << endl << grid.size() << endl;
}
void bsplines::printvals(uint m, uint n)
{
cout << "content in spline (B" << m << "," << n << ") vals vector: " << endl;
for (auto& el : splines[n][m].vals)
cout << el << " ";
cout << endl;
}
void bsplines::writesplines()
{
for (uint o = 0; o < order; o++)
for (auto& sp : splines[o])
{
uint i = sp.index;
ostringstream namestream;
namestream << "B(" << fixed << setprecision(1) << i << ","
<< fixed << setprecision(1) << o << ").csv";
string filename = namestream.str();
ofstream fs;
fs.open(filename);
if (o < order - 1)
{
for (uint k = 0; k < sp.vals.size(); k++)
fs << sp.vals[k] << "," << 0 << "," << 0 << endl;
fs.close();
}
else
{
for (uint k = 0; k < sp.vals.size(); k++)
fs << sp.vals[k] << "," << sp.d1[k] << "," << sp.d2[k] << endl;
fs.close();
}
cout << "write " << sp.vals.size() << " numbers to " << filename << endl;
}
}
编辑：已更新。h-file。
解决了。你们非常亲密。您所要做的就是扩展衍生品 与
1 / (x_max - x_min)（一阶导数）1 / (x_max - x_min)^2（二阶导数）。
TLDR：您在拟合样条曲线时将 x 值归一化为介于 0 和 1 之间，但没有缩放 y 值。
您实际拟合的不是样条拟合 x^2，而是：
x_norm = (x - x_min) / (x_max - x_min)
y = x_norm**2
因此使用链式法则 y = x_norm**2 的一阶导数为 2x / (x_max - x_min)，二阶导数为 2 / (x_max - x_min)**2.
完整示例代码：
#include <iostream>
#include <Eigen/Core>
#include <unsupported/Eigen/Splines>
using namespace Eigen;
using namespace std;
VectorXd normalize(const VectorXd &x) {
VectorXd x_norm;
x_norm.resize(x.size());
const double min = x.minCoeff();
const double max = x.maxCoeff();
for (int k = 0; k < x.size(); k++) {
x_norm(k) = (x(k) - min)/(max - min);
}
return x_norm;
}
int main() {
typedef Spline<double, 1, 3> Spline1D;
typedef SplineFitting<Spline1D> Spline1DFitting;
const Vector4d x{0, 1, 2, 4};
const Vector4d y = (x.array().square()); // x^2
const auto knots = normalize(x); // Normalize x to be between 0 and 1
const double scale = 1 / (x.maxCoeff() - x.minCoeff());
const double scale_sq = scale * scale;
Spline1D spline = Spline1DFitting::Interpolate(y.transpose(), 3, knots);
cout << "1st deriv at x = 0: " << spline.derivatives(0.00, 1)(1) * scale << endl;
cout << "1st deriv at x = 1: " << spline.derivatives(0.25, 1)(1) * scale << endl;
cout << "1st deriv at x = 2: " << spline.derivatives(0.50, 1)(1) * scale << endl;
cout << "1st deriv at x = 3: " << spline.derivatives(0.75, 1)(1) * scale << endl;
cout << "1st deriv at x = 4: " << spline.derivatives(1.00, 1)(1) * scale << endl;
cout << endl;
/**
* IMPORTANT NOTE: Eigen's spline module is not documented well. Once you fit a spline
* to find the derivative of the fitted spline at any point u [0, 1] you call:
*
* spline.derivatives(u, 1)(1)
* ^ ^ ^
* | | |
* | | +------- Access the result
* | +---------- Derivative order
* +------------- Parameter u [0, 1]
*
* The last bit `(1)` is if the spline is 1D. And value of `1` for the first
* order. `2` for the second order. Do not forget to scale the result.
*
* For higher dimensions, treat the return as a matrix and grab the 1st or
* 2nd column for the first and second derivative.
*/
cout << "2nd deriv at x = 0: " << spline.derivatives(0.00, 2)(2) * scale_sq << endl;
cout << "2nd deriv at x = 1: " << spline.derivatives(0.25, 2)(2) * scale_sq << endl;
cout << "2nd deriv at x = 2: " << spline.derivatives(0.50, 2)(2) * scale_sq << endl;
cout << "2nd deriv at x = 3: " << spline.derivatives(0.75, 2)(2) * scale_sq << endl;
cout << "2nd deriv at x = 4: " << spline.derivatives(1.00, 2)(2) * scale_sq << endl;
return 0;
}
示例输出：
1st deriv at x = 0: 4.52754e-16
1st deriv at x = 1: 2
1st deriv at x = 2: 4
1st deriv at x = 3: 6
1st deriv at x = 4: 8
2nd deriv at x = 0: 2
2nd deriv at x = 1: 2
2nd deriv at x = 2: 2
2nd deriv at x = 3: 2
2nd deriv at x = 4: 2