Commit 29e6095b authored by Nils Rummler's avatar Nils Rummler
Browse files

solution of bandmatrix is incorrect

parent 92ad1b57
......@@ -41,7 +41,7 @@ namespace ITABase
public:
CPiecewisePolynomial( const std::vector<float>& vdBreakPoints, const std::vector<float>& vdCoeffs, const int iOrder );
// Returns the polynomial order
//! Returns the polynomial order
inline int Order() const { return iOrder; };
//! Returns the number of intervals (pieces) between break points
inline int NumIntervals() const { return vdBreakPoints.size() - 1; };
......@@ -64,7 +64,7 @@ namespace ITABase
std::vector<float>::const_iterator CoefficientsOfInterval(const int idxInterval) const;
private:
// Order of polynomials
//! Order of polynomials
int iOrder;
//! Vector with break points. Defines the interval for each polynomial.
std::vector<float> vdBreakPoints;
......
......@@ -21,6 +21,7 @@ CPiecewisePolynomial::CPiecewisePolynomial(const std::vector<float>& vdBreakPoin
const int nCoeffs = NumIntervals() * (iOrder + 1);
if (vdCoeffs.size() != nCoeffs)
ITA_EXCEPT_INVALID_PARAMETER("Invalid number of coefficients: Expected " + std::to_string(nCoeffs) + " but got " + std::to_string(vdCoeffs.size()) + ".");
}
float CPiecewisePolynomial::Value(const float xValue) const
......@@ -72,6 +73,87 @@ std::vector<float>::const_iterator ITABase::CPiecewisePolynomial::CoefficientsOf
CPiecewisePolynomial ITABase::Spline(const std::vector<float>& vdSupportingPoints, const std::vector<float>& vdDataPoints, const int iPolynomialOrder)
{
//TODO: Calculate coefficients for polynomials
return CPiecewisePolynomial(vdSupportingPoints, { 0.0, 0.0 }, iPolynomialOrder);
/*setUp the bandmatrix with the the h values, called M. The vektor with the Koefficents c is called c.
The vector with the linear approximations on the right hand side is called g.
M * c = g
Solve the bandmatrix M with the thomsas algorithm and reconstruct the other coefficients according to:
www.tm-mathe.de/Themen/html/funnatsplines.html
out of the c vector.
*/
//Calculate h_i
int iPoints = vdSupportingPoints.size();
std::vector<float> vdH(iPoints-1);
for (int i = 0; i < iPoints-1; i++) {
vdH[i] = vdSupportingPoints[i + 1] - vdSupportingPoints[i];
}
//SetUp the band matrix with lower, middle and upper diagonals. momentary only natural splines
std::vector<float> vdLower(iPoints - 3), vdMiddle(iPoints - 2), vdUpper(iPoints - 3);
for (int i = 0; i < iPoints - 3;i++) {
vdLower[i] = vdH[i + 1];
}
vdUpper = vdLower;//TODO: non natural splines behave different -> put it in the loop
for (int i = 0; i < iPoints - 2;i++) {
vdMiddle[i] = 2 * (vdH[i] + vdH[i + 1]);
}
//excitation vector
//TODO only for natural splines with S0=S_iPoints = 0
std::vector<float> vdG(iPoints - 2);
for (int i = 0; i < iPoints - 2; i++) {
vdG[i] = (vdDataPoints[i + 2] - vdDataPoints[i + 1]) / vdH[i + 1] - (vdDataPoints[i + 1] - vdDataPoints[i]) / vdH[i];
vdG[i] *= 3;
}
// solving the band matrix with Thomas Algorithm
//https://www.quantstart.com/articles/Tridiagonal-Matrix-Solver-via-Thomas-Algorithm/
std::vector<float> vdC_star(iPoints - 3);
std::vector<float> vdG_star(iPoints - 2);
vdC_star[0] = vdUpper[0] / vdMiddle[0];
vdG_star[0] = vdG[0] / vdMiddle[0];
for (int i = 1; i < iPoints - 2; i++) {
if(i<iPoints-3){
vdC_star[i] = vdUpper[i] / (vdMiddle[i] - vdLower[i - 1] * vdC_star[i - 1]);
}
vdG_star[i] = (vdG[i] - vdG_star[i - 1] * vdLower[i - 1]) / (vdMiddle[i] - vdC_star[i - 1] * vdLower[i - 1]);
}
std::vector<float> vdB(iPoints - 1); // coeffzent b and solution of bandmatrix
vdB[iPoints - 2] = vdG_star[iPoints - 3];
for (int i = iPoints - 3; i > 0;) {
vdB[i] = vdG_star[i] - vdC_star[i-1] * vdB[i + 1];
i--;
}
//bandmatrix is solved, now get the other coefficients
std::vector<float> vdA(iPoints - 1);
std::vector<float> vdC(iPoints - 1);
std::vector<float> vdD(iPoints - 1);
for (int i = 0; i < iPoints - 2; i++) {
vdC[i] = (vdDataPoints[i + 1] - vdDataPoints[i]) / vdH[i] - 1 / 3 * vdH[i] * (2 * vdB[i] + vdB[i + 1]);
vdA[i] = (vdB[i + 1] - vdB[i]) / (3 * vdH[i]);
}
vdD = vdDataPoints;
vdD.pop_back();
int iCoefficients = (iPoints - 1) * (iPolynomialOrder + 1);
std::vector<float> vdCoefficients(iCoefficients);
int k = 0;
for (int i = 0; i < iCoefficients;) {
vdCoefficients[i] = vdA[k];
vdCoefficients[i + 1] = vdB[k];
vdCoefficients[i + 2] = vdC[k];
vdCoefficients[i + 3] = vdD[k];
i += 4;
k += 1;
}
return CPiecewisePolynomial(vdSupportingPoints, vdCoefficients, iPolynomialOrder);
}
......@@ -21,11 +21,21 @@
#include <ITABase/Spline.h>
//#include <iostream>
#include <iostream>
using namespace std;
//using namespace std;
int main(int iNumInArgs, char* pcInArgs[])
{
std::vector<float> vdSupportingPoints{1,7,12,15,19};
std::vector<float> vdDataPoints = {8,10,7,8,7};
std::vector<float> vdSolution{-0.01045, 0, 0.7091, 8, 0.0304, -1.88185, -0.4194, 10, -0.05005, 0.268014, -0.02026, 7, 0.01520, -0.182432, 0.2365, 8}
int iPolynomialOrder = 3;
IW_ITA_BASE_SPLINE::ITABase::CPiecewisePolynomial myPoly = IW_ITA_BASE_SPLINE::ITABase::Spline(vdSupportingPoints,vdDataPoints,iPolynomialOrder);
/*for (int i=0; it < vdDataPoints.size(); i++) {
std::cout << vdSolution[i] - myPoly.get_Coefficient(i);
}
*/
return 0;
}
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