Commit 4ef7d6b0 authored by Nils Rummler's avatar Nils Rummler
Browse files

still only third order, but with possibility to get values from approximation

parent e67dcbd7
......@@ -52,6 +52,8 @@ namespace ITABase
//! Evaluate the piecewise polynomial at a given x value.
float Value(const float xValue) const;
//! Finds the index of Interval at a given x value
int IntervalIndex(const float xValue) const;
//! Return a piecewise polynomial containing the derivations of the polynomials of this one
CPiecewisePolynomial Derivation() const;
......
......@@ -3,6 +3,7 @@
#include <string>
#include <algorithm>
#include <cmath>
using namespace ITABase;
......@@ -26,18 +27,27 @@ CPiecewisePolynomial::CPiecewisePolynomial(const std::vector<float>& vdBreakPoin
float CPiecewisePolynomial::Value(const float xValue) const
{
//TODO: Check if x-value is inside boundaries
//TODO: Find index of interval, get corresponding coefficients and calculate values
// Check if x-value is inside boundaries
// Find index of interval
int iInterval = IntervalIndex(xValue);
//TODO: get corresponding coefficients and calculate values
std::vector<float>::const_iterator itCoef = CoefficientsOfInterval(iInterval);
//f(x) = a1(x−x1)^3 + b1(x−x1)^2 + c1(x−x1) + d1
float x0 = vdBreakPoints[iInterval];
if (iOrder == 3)
return itCoef[0] * pow(xValue - x0, 3) + itCoef[1] * pow(xValue - x0, 2) + itCoef[2] * (xValue - x0) + itCoef[3];
else
ITA_EXCEPT_INVALID_PARAMETER("currently only splines of order 3");
}
//Something like this might be necessary
std::vector<float>::const_iterator itStart = CoefficientsOfInterval(0);
std::vector<float>::const_iterator itEnd = itStart + NumCoefficients();
for (std::vector<float>::const_iterator it = itStart; itStart < itEnd; it++)
{
int ITABase::CPiecewisePolynomial::IntervalIndex(const float xValue) const
{
if (xValue < vdBreakPoints[0] || xValue > vdBreakPoints[NumIntervals()])
ITA_EXCEPT_INVALID_PARAMETER("xValue is out of bounds");
for (int i = 1; i < vdBreakPoints.size(); i++) {
if (vdBreakPoints[i] > xValue)
return i - 1;
}
return 0.0;
}
CPiecewisePolynomial CPiecewisePolynomial::Derivation() const
......@@ -80,12 +90,15 @@ CPiecewisePolynomial ITABase::Spline(const std::vector<float>& vdSupportingPoint
https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&cad=rja&uact=8&ved=2ahUKEwjcofWO_IrqAhUux4sKHT52AigQFjABegQIBRAB&url=http%3A%2F
out of the b vector.
*/
if (vdSupportingPoints.size()!= vdDataPoints.size())
ITA_EXCEPT_INVALID_PARAMETER("vdSupportingPoints and vdDataPoints need to have the same size");
//TODO andere als Ordnung 3
//Calculate h_i
//Calculate h_i, formular 1.27b
int iPoints = vdSupportingPoints.size();
std::vector<float> vdH(iPoints-1);
for (int i = 0; i < iPoints-1; i++) {//TODO needs tp be positiv
for (int i = 0; i < iPoints-1; i++) {
vdH[i] = vdSupportingPoints[i + 1] - vdSupportingPoints[i];
if (vdH[i] < 0) {
ITA_EXCEPT_INVALID_PARAMETER("vdSupportingPoints needs to be sorted in ascending order");
......@@ -106,11 +119,10 @@ CPiecewisePolynomial ITABase::Spline(const std::vector<float>& vdSupportingPoint
//excitation vector
//TODO only for natural splines with S0=S_iPoints = 0
std::vector<float> vdR(iPoints - 2);
for (int i = 0; i < iPoints - 2; i++) {
vdR[i] = (vdDataPoints[i + 2] - vdDataPoints[i + 1]) / vdH[i + 1] - (vdDataPoints[i + 1] - vdDataPoints[i]) / vdH[i];
vdR[i] *= 3;
}
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.");
......@@ -148,13 +160,11 @@ CPiecewisePolynomial ITABase::Spline(const std::vector<float>& vdSupportingPoint
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;
for (int k = 0; k < iPolynomialOrder + 1;k++) {
vdCoefficients[4*k] = vdA[k];
vdCoefficients[4*k + 1] = vdB[k];
vdCoefficients[4*k + 2] = vdC[k];
vdCoefficients[4*k + 3] = vdD[k];
}
return CPiecewisePolynomial(vdSupportingPoints, vdCoefficients, iPolynomialOrder);
......
......@@ -32,7 +32,7 @@ int main(int iNumInArgs, char* pcInArgs[])
//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);
std::cout << "Quick and dirty test: "<< myPoly.Value(5)<<" -hier soll grob 10 rauskommen;)";
/*for (int i=0; it < vdDataPoints.size(); i++) {
std::cout << vdSolution[i] - myPoly.get_Coefficient(i);
......
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