K
So, hier mal meine Splineklassen, die wahrscheinlich nicht höchsten Programmierprofiansprüchen genügen, aber dennoch funktionieren...bei mir zumindest :).
Folgende Klassen/Includes werden verwendet, die ich jetzt nicht extra posten möchte, da das ganze sonst zu unübersichtlich wird.
CVektor: Eine Vektorklasse von mir abgeleitet von std::vector. Die mathematischen Operatoren sind überladen.
CMatrix: Meine Matrixklasse zum Lösen von Gleichungssystemen etc.
CFehler: Meine Fehlerklasse.
allgemein.h: Eine Includedatei, in der alle benötigten Dateien aus der STL eingebunden werden.
So, jetzt zu den eigentlichen Klassen.
Zuerst die Klasse CSpline (repräsentiert ein Polynom):
#ifndef _spline_h
#define _spline_h
#include "vektor.h"
using namespace std;
class CSpline
{
private:
//Private Variablen
//Faktoren
CVektor<double> ffaktoren;
//Grad des Polynoms
const unsigned fgrad;
//Grenzen des Bereiches, in dem das Spline gelten soll
pair<double, double> fbereich;
public:
//Konstruktoren
CSpline();
CSpline(const double & axvon, const double & axbis, const unsigned & agrad);
//Copy-Konstruktor
CSpline(const CSpline &);
//Destruktor
virtual ~CSpline();
const double gib_funktionswert(const double & ax) const;
const CVektor<double> gib_funktionswerte(const CVektor<double> & axvec) const;
void setze_faktor(const unsigned & anr, const double & awert) { ffaktoren.at(anr) = awert; };
void setze_faktoren(const CVektor<double> & awerte) { ffaktoren = awerte; };
const CVektor<double> & gib_faktoren() const { return ffaktoren; };
const unsigned & gib_grad() const { return fgrad; };
const pair<double, double> & gib_bereich() const { return fbereich; };
void setze_bereich(const pair<double, double> & abereich) { fbereich = abereich; };
};
#endif
#include "allgemein.h"
#include "spline.h"
#include "vektor.h"
CSpline::CSpline():fgrad(3), fbereich(make_pair(.0, .0))
{
throw CFehler("Standardkonstruktor bnicht verwenden!","CSpline()","CSpline");
}
CSpline::CSpline(const double & axvon, const double & axbis, const unsigned & agrad):fgrad(agrad), fbereich(make_pair(axvon, axbis))
{
ffaktoren.resize(agrad+1);
}
CSpline::CSpline(const CSpline & aspline):fgrad(aspline.gib_grad()), fbereich(aspline.gib_bereich())
{
ffaktoren=aspline.gib_faktoren();
}
CSpline::~CSpline()
{
ffaktoren.clear();
}
const double CSpline::gib_funktionswert(const double & ax) const
{
double erg = ffaktoren.at(0);
for(int i=1; i<=(int)fgrad; i++)
{
erg+=ffaktoren.at(i)*pow(ax, (int)i);
}
return erg;
}
const CVektor<double> CSpline::gib_funktionswerte(const CVektor<double> & axvec) const
{
CVektor<double> erg;
for(unsigned i=0; i< axvec.size(); i++)
{
erg.push_back(gib_funktionswert(axvec.at(i)));
}
return erg;
}
Nun die Klasse CSplineApprox, die sich um das Lösen des linearen Gleichungssystems kümmert:
#ifndef _splineapprox_h
#define _splineapprox_h
#include "vektor.h"
#include "matrix.h"
class CSpline;
using namespace std;
class CSplineApprox
{
private:
//Private Variablen
//Soll beim ersten und letzten Punkt die erste oder zweite Ableitung = 0 sein?
//2: zweite Ableitung
//1: erste Ableitung
int fgrenzbedingung;
const int fsplinegrad;
//Alle einzelnen Splines
vector<CSpline*> fsplines;
//Stuetzstellen auf der x-Achse
CVektor<double> fx;
//Stuetzfunktionswerte
CVektor<double> fy;
public:
//Konstruktoren
CSplineApprox(const int & agrenzbedingung = 2);
//Copy-Konstruktor
CSplineApprox(const CSplineApprox &);
//Destruktor
virtual ~CSplineApprox();
//void setze_x_stuetzstellen(const CVektor<double> & ax) { fx = ax; };
const CVektor<double> & gib_x_stuetzstellen() const { return fx; };
//void setze_y_stuetzfunktionswerte(const CVektor<double> & ay) { fy = ay; };
const CVektor<double> & gib_y_stuetzfunktionswerte() const { return fy; };
const vector<CSpline*> & gib_splines() const { return fsplines; };
//Approximierte Funktion fuer vorgegebene x-Werte berechnen
const CVektor<double> gib_funktionswerte(const CVektor<double> & ax);
//Splines anhand der Anzahl der Stuetzstellen initialisieren
void init_splines(const CVektor<double> & ax, const CVektor<double> & ay);
void approximiere();
void setze_grenzbedingung(const int & awert) { fgrenzbedingung = awert; };
const int & gib_grenzbedingung() const { return fgrenzbedingung; };
};
#endif
#include "allgemein.h"
#include "splineapprox.h"
#include "vektor.h"
#include "spline.h"
CSplineApprox::CSplineApprox(const int & agrenzbedingung):fsplinegrad(3), fgrenzbedingung(agrenzbedingung)
{
}
CSplineApprox::CSplineApprox(const CSplineApprox & aspline):fsplinegrad(3)
{
}
CSplineApprox::~CSplineApprox()
{
}
void CSplineApprox::init_splines(const CVektor<double> & ax, const CVektor<double> & ay)
{
fx=ax;
fy=ay;
//Splines anlegen
for(int i=1; i<(int)fx.size(); i++)
{
//Kubische Splines anlegen
fsplines.push_back(new CSpline(ax.at(i-1), ax.at(i), fsplinegrad));
}
}
void CSplineApprox::approximiere()
{
const int & anz_faktoren = fsplinegrad+1;
const int & matrixdim = ((int)fx.size()-1)*anz_faktoren;
const int & halbematrixdim = matrixdim/2;
//Matrix und Vektor der rechten Seite fuellen
CMatrix<double> A(matrixdim, matrixdim);
CVektor<double> b(matrixdim);
//Matrix mit Nullen fuellen
A.init(.0);
//Schleife ueber alle Zeilen
for(int izeile=0; izeile<matrixdim; izeile++)
{
//Vektor der rechten Seite setzen
if(izeile<halbematrixdim)
{
b[izeile]=fy[(izeile+1)/2];
}
else
{
b[izeile]=.0;
}
if(izeile<matrixdim/2)
{
//Gleichungen der Splinefunktion
const int & von = (izeile/2)*anz_faktoren;
const int & bis = von+anz_faktoren-1;
for(int ispalte=von;ispalte<=bis; ispalte++)
{
const double & wert = pow(fx[(izeile+1)/2], ispalte%anz_faktoren);
A.setze_element(izeile, ispalte, wert);
}
}
else if(izeile>=halbematrixdim && izeile<halbematrixdim+(halbematrixdim-2)/2)
{
//Gleichungen der 1. Ableitung der Splinefunktion
const int & von = (izeile-halbematrixdim)*anz_faktoren+1;
const int & bis = von+anz_faktoren-2;
for(int ispalte=von;ispalte<=bis; ispalte++)
{
const int & faktor = ispalte-anz_faktoren*(izeile-halbematrixdim);
const double & wert = faktor*pow(fx[izeile-halbematrixdim+1], ispalte%anz_faktoren-1);
A.setze_element(izeile, ispalte, wert);
A.setze_element(izeile, ispalte+anz_faktoren, -wert);
}
}
else if(izeile>=halbematrixdim+(halbematrixdim-2)/2 && izeile<matrixdim-2)
{
//Gleichungen der 2. Ableitung der Splinefunktion
const int & von = (izeile-halbematrixdim-(halbematrixdim-2)/2)*anz_faktoren+2;
const int & bis = von+anz_faktoren-3;
for(int ispalte=von;ispalte<=bis; ispalte++)
{
const int & faktor = fak(ispalte-anz_faktoren*(izeile-halbematrixdim-(halbematrixdim-2)/2));
const double & wert = faktor*pow(fx[izeile-3*halbematrixdim/2+2], ispalte%anz_faktoren-2);
A.setze_element(izeile, ispalte, wert);
A.setze_element(izeile, ispalte+anz_faktoren, -wert);
}
}
else if(izeile==matrixdim-2)
{
if(fgrenzbedingung==2)
{
//Gleichungen der 2. Ableitung des ersten und letzten Stuetzwertes
const int & von = 2;
const int & bis = von+anz_faktoren-3;
for(int ispalte=von;ispalte<=bis; ispalte++)
{
const int & faktor = fak(ispalte);
const double & wert1 = faktor*pow(fx[0], ispalte%anz_faktoren-2);
const double & wert2 = faktor*pow(*(fx.end()-1), ispalte%anz_faktoren-2);
A.setze_element(izeile, ispalte, wert1);
A.setze_element(izeile+1, matrixdim-anz_faktoren+ispalte, wert2);
}
}
else if(fgrenzbedingung==1)
{
//Gleichungen der 1. Ableitung des ersten und letzten Stuetzwertes
const int & von = 1;
const int & bis = von+anz_faktoren-2;
for(int ispalte=von;ispalte<=bis; ispalte++)
{
const int & faktor = ispalte;
const double & wert1 = faktor*pow(fx[0], ispalte%anz_faktoren-1);
const double & wert2 = faktor*pow(*(fx.end()-1), ispalte%anz_faktoren-1);
A.setze_element(izeile, ispalte, wert1);
A.setze_element(izeile+1, matrixdim-anz_faktoren+ispalte, wert2);
}
}
else
throw CFehler("Wert fuer Grenzbedingung nicht gueltig.", "approximiere", "CSplineApprox");
}
}
//Gleichungssystem loesen
CVektor<double> x = loese_lgs(A, b);
//Faktoren der einzelnen Splines setzen
for(int i=0; i<(int)x.size(); i++)
{
fsplines[i/(fsplinegrad+1)]->setze_faktor(i%(fsplinegrad+1), x[i]);
}
}
const CVektor<double> CSplineApprox::gib_funktionswerte(const CVektor<double> & ax)
{
CVektor<double> erg;
vector<CSpline*>::const_iterator spline_it = fsplines.begin();
double splinegrenze = (*spline_it)->gib_bereich().second;
for(unsigned i=0; i<ax.size(); i++)
{
//richtiges Spline finden
while(splinegrenze<ax[i])
{
spline_it++;
splinegrenze = (*spline_it)->gib_bereich().second;
}
//Funktionswert speichern
erg.push_back((*spline_it)->gib_funktionswert(ax[i]));
}
return erg;
}
So, das war´s auch schon...