So erstellen Sie eine rekursive Spline-Funktion in C ++

Im Moment arbeite ich an einer Differentialgleichungslösungsmethode namens Basis-Spline-Kollokation. Ich habe Probleme damit, eine Methode zum Erstellen eines Splines beliebiger Ordnung mit der Beziehung zu erstellen. mit dem Anfangsbedingung und ich habe Probleme, selbst mit diesem Problem zu beginnen, da es rekursiv ist und entweder von "oben" oder "unten" beginnen kann, und ich stoße auf einen allgemeinen Schreibblocktyp von der Sache, wo ich mich nicht auf das konzentrieren kann, was ich tun muss.B 1 i (x)={

$$B^{k+1}_{i}(x)= \frac{x-x_i}{x_{k+i}-x_i}B^k_i + \frac{x_{k+i+1}-x}{x_{k+i+1}-x_{i+1}}B^k_{i+1}(x)$$ $$B^1_i(x)=\begin{cases} 1 & \; \text{for } \; x_i \leq x < x_{i+1} \\ 0 & \; \text{otherwise} \end{cases}$$

Ich kann empfehlen, das NURBS-Buch zu konsultieren , das ein klassischer Text zu diesem Thema zu sein scheint. Der Algorithmus selbst ist auf Seite 72 angegeben und kann online angezeigt werden.

Ich unterstütze das NURBS-Buch und möchte auch hervorheben, dass diese rekursive Formel zwar normalerweise die Art und Weise darstellt, wie man die B-Spline-Basis ausdrückt (wie in einem Artikel), aber nicht, wie Sie die Basisfunktionen implementieren. Da alle Funktionen höherer Ordnung auf niedrigerer Ordnung basieren, können Sie alle Funktionen ungleich Null gleichzeitig berechnen und die Auswertungen niedrigerer Ordnung wiederverwenden. So macht es der Piegl-Algorithmus.$p+1$

Ich weiß ehrlich gesagt nicht, wie effizient dies ist, aber eine Möglichkeit, dies zu tun, sind C ++ - Vorlagen:

Die Reihenfolge ist k, t ist die Knotenstruktur und x ist der gewünschte Wert.

template <int k> 
real BSpline(real x, real *t)
{
    if (*t <= x && x < *(t+k))
    {
        real a = (x - *t) / (*(t+k-1) - *t);
        real b = (*(t+k) - x) / (*(t+k) - *(t+1));

        return a * BSpline<k-1>(x, t) + b * BSpline<k-1>(x, (t+1));
    }
    else
        return 0;
};

template <>
real BSpline<1>(real x, real *t)
{
    if (*t <= x && x < *(t+1))
        return 1.;
    else
        return 0.;
};