Commit cfcc03d0 authored by Nils Rummler's avatar Nils Rummler
Browse files

Bandmatrix solver is now an independent function within the spline header

parent ee88af43
......@@ -87,6 +87,9 @@ namespace ITABase
//! Uses cubic splines to create a piecewise polynomial (order 3) for the given data pairs. vdSupportingPoints need to be sorted in ascending order
ITA_BASE_API CPiecewisePolynomial CubicSpline(const std::vector<float>& vdSupportingPoints, const std::vector<float>& vdDataPoints) { return Spline(vdSupportingPoints, vdDataPoints, 3); };
//! Solves a Tridiagonal Bandmatrix with Thomas algorithm. Expects the 3 diagonals (lower,middel and upper) and the excitation vector.
std::vector<float> BandmatrixSolver(const std::vector<float> vdLower,const std::vector<float> vdMiddle,const std::vector<float> vdUpper,const std::vector<float> vdExcitation);
}
#endif // IW_ITA_BASE_SPLINE
......@@ -130,26 +130,9 @@ CPiecewisePolynomial ITABase::Spline(const std::vector<float>& vdSupportingPoint
for (int i = 0; i < iPoints - 2; i++)
vdR[i] = 3 *( (vdDataPoints[i + 2] - vdDataPoints[i + 1]) / vdH[i + 1] - (vdDataPoints[i + 1] - vdDataPoints[i]) / vdH[i]);
//TODO auslagern
//solving bandmatrix with "numerical recipees in C" page 51
//if (vdMiddle[0] == 0)
//ITA_EXCEPT_INVALID_PARAMETER("Matrix needs to be semi positiv definit.");
std::vector<float> vdGamma(iPoints - 2);
float fBeta = vdMiddle[0];
std::vector<float> vdB(iPoints - 1); // coeffzent b and solution of bandmatrix
vdB[1] = vdR[0]/vdMiddle[0];
for (int j = 1; j < iPoints - 2; j++) {
vdGamma[j] = vdUpper[j - 1] / fBeta;
fBeta = vdMiddle[j] - vdLower[j - 1] * vdGamma[j];
vdB[j + 1] = (vdR[j] - vdLower[j - 1] * vdB[j]) / fBeta;
}
for (int j = iPoints - 3; j > 0; j--) {//backsubstitution
vdB[j] -= vdGamma[j] * vdB[j + 1];
}
//bandmatrix is solved, now get the other coefficients
//solve bandmatrix and get the other coefficients
std::vector<float> vdB = BandmatrixSolver(vdLower,vdMiddle,vdUpper,vdR);
vdB.insert(vdB.begin(), 0.0);
std::vector<float> vdA(iPoints - 1);
std::vector<float> vdC(iPoints - 1);
......@@ -177,3 +160,31 @@ CPiecewisePolynomial ITABase::Spline(const std::vector<float>& vdSupportingPoint
return CPiecewisePolynomial(vdSupportingPoints, vdCoefficients, iPolynomialOrder);
}
std::vector<float> ITABase::BandmatrixSolver(const std::vector<float> vdLower,const std::vector<float> vdMiddle, const std::vector<float> vdUpper, const std::vector<float> vdExcitation)
{
//solving bandmatrix with "numerical recipees in C" page 51
if (vdMiddle[0] == 0)
ITA_EXCEPT_INVALID_PARAMETER("Matrix needs to be semi positiv definit.");
if (vdMiddle.size() != vdExcitation.size() || vdLower.size() != vdUpper.size() || vdExcitation.size() != vdLower.size() + 1)
ITA_EXCEPT_INVALID_PARAMETER("Dimensions of arguments missmatch.");
int n = vdMiddle.size();
std::vector<float> vdGamma(n);
float fBeta = vdMiddle[0];
std::vector<float> vdB(n); //solution of bandmatrix
vdB[0] = vdExcitation[0] / vdMiddle[0];
for (int j = 1; j < n; j++) {
vdGamma[j] = vdUpper[j - 1] / fBeta;
fBeta = vdMiddle[j] - vdLower[j - 1] * vdGamma[j];
vdB[j] = (vdExcitation[j] - vdLower[j - 1] * vdB[j-1]) / fBeta;
}
for (int j = n-2; j >= 0; j--) {//backsubstitution
vdB[j] -= vdGamma[j+1] * vdB[j + 1];
}
return vdB;
}
......@@ -23,10 +23,11 @@ Coeffs = [-24.4868088000000 0 27.4868088000000 77;...%Copy & Paste from the resu
-29.6697369000000 89.0092087000000 -34.3394737000000 26];
myPoly = mkpp(X,Coeffs);
%% Plot
xx = [-10:0.1:10];
plot(X,Y,'o',xx,spline(X,Y,xx));
hold on
grid on
title 'Splines: Matlab vs. ITA'
plot(xx,ppval(myPoly2,xx));
plot(xx,spline(X,Y,xx)-ppval(myPoly2,xx));
plot(xx,ppval(myPoly,xx));
plot(xx,spline(X,Y,xx)-ppval(myPoly,xx));
legend('Data','Matlab','ITA','Diff');
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment