Natural Cubic Spline

Im Artikel Gaußverfahren mit Pivoting hatte ich bereits darauf hingewiesen, dass mit diesem Verfahren auch Splines berechnet werden können. Hier ist die VBA-Funktion zum Erstellen von natürlichen kubischen Splines (natural cubic splines).


Der Algorithmus wurde in Anlehnung an den Wikipedia Artikel zum Thema Splines umgesetzt. Der Aufruf erfolgt, zum Beispiel, so:

VB-Code:
Sub TestSpline()
    Dim varX As Variant, varXX As Variant, varY As Variant
    Dim dblX() As Double, dblXX() As Double, dblY() As Double
    Dim i As Long

    With Sheets("NaturalCubicSpline")
        varX = .Range("A2:A9")
        varY = .Range("B2:B9")
        varXX = .Range("D2:D22")
        ReDim dblX(UBound(varX), 1)
        ReDim dblY(UBound(varX), 1)
        ReDim dblXX(UBound(varXX), 1)
        For i = 1 To UBound(varX)
            dblX(i, 1) = varX(i, 1)
            dblY(i, 1) = varY(i, 1)
        Next

        For i = 1 To UBound(varXX)
            dblXX(i, 1) = varXX(i, 1)
        Next

        .Range("E2:E22") = NaturalCubicSpline(dblX, dblY, dblXX)
    End With
End Sub

Code eingefügt mit Syntaxhighlighter 4.15

In der aktuellen Version werden, wie bereits beim Algorithmus zum Gauß-Verfahren, immer zweidimensionale Arrays vom Datentyp double benötigt. Vielleicht werde ich alles in nächster Zeit auch Variant umstellen.

Quellcode

Der Quellcode, siehe unten, kann auch als
Modul heruntergeladen werden.

Option Explicit
Option Base 1

'------------------------------------------------------------------
' Module    : NaturalCubicSpline
' DateTime  : 10.09.2008 13:23
' Author    : Tobias - tobiasschmid.de - vba-blog.de
' Purpose   : Calculate an natural cubic spline
'------------------------------------------------------------------

'------------------------------------------------------------------
' Procedure : TestSplineMitVariant
' DateTime  : 10.09.2008 13:24
' Author    : Tobias - tobiasschmid.de - vba-blog.de
' Purpose   : Example of usage
'------------------------------------------------------------------
Sub TestSpline()
    Dim varX As Variant, varXX As Variant, varY As Variant

    With Sheets("NaturalCubicSpline")
        varX = .Range("A2:A9")
        varY = .Range("B2:B9")
        varXX = .Range("D2:D22")
        .Range("E2:E22") = NaturalCubicSpline(varX, varY, varXX)
    End With
End Sub

'------------------------------------------------------------------
' Procedure : NaturalCubicSpline
' DateTime  : 10.09.2008 13:24
' Author    : Tobias - tobiasschmid.de - vba-blog.de
' Purpose   : Calculating an natural cubic spline
'------------------------------------------------------------------
Public Function NaturalCubicSpline(x As Variant, y As Variant, xx As VariantAs Variant
    Dim Matrix() As Variant, Vektor() As Variant, z() As Double, yy() As Double, h() As Double
    Dim i As Long, j As Long
    ReDim Matrix(UBound(x), UBound(x))
    ReDim Vektor(UBound(x), 1)
    ReDim z(UBound(x), 1)
    ReDim h(UBound(x), 1)
    ReDim yy(UBound(xx), 1)

    For i = 1 To UBound(x) - 1
        h(i, 1) = x(i + 1, 1) - x(i, 1)
    Next
    'z_0 = 0
    'h_i-1 z_i-1 + 2(h_i-1 + h_i) z_i + h_i z_i+1 = 6((y_i+1 -y_i)/h_i - (y_i -y_i-1)/h_i-1)
    'z_n = 0

    Matrix(1, 1) = 1
    For i = 2 To UBound(x) - 1
        Vektor(i, 1) = 6 * ((y(i + 1, 1) - y(i, 1)) / h(i, 1) - (y(i, 1) - y(i - 1, 1)) / h(i - 1, 1))
        Matrix(i, i - 1) = h(i - 1, 1)
        Matrix(i, i) = 2 * (h(i - 1, 1) + h(i, 1))
        Matrix(i, i + 1) = h(i, 1)
    Next
    Matrix(UBound(x), UBound(x)) = 1

    For i = 1 To UBound(Matrix)
        For j = 1 To UBound(Matrix)
            'Debug.Print Matrix(i, j),
        Next
        'Debug.Print ""
    Next

    z = GaussEliminationWithPivoting(Matrix, Vektor)
    For i = 1 To UBound(Matrix)
        'Debug.Print z(i, 1)
    Next
    'S_i(x) = (z_i+1(x-x_i)^3+z_i(x_i+1-x)^3)/(6 h_i)+ _
     (y_i+1/h_i - h_i/6 z_i+1)(x-x_i) + _
     (y_i/h_i - h_i/6 z_i) (x_i+1-x)

    j = 1
    For i = 1 To UBound(xx)
        'Do While Not (xx(i, 1) <= x(j + 1, 1))
        Do While xx(i, 1) > x(j + 1, 1)
            j = j + 1
        Loop
        yy(i, 1) = (z(j + 1, 1) * (xx(i, 1) - x(j, 1)) ^ 3 + z(j, 1) * (x(j + 1, 1) - xx(i, 1)) ^ 3) / (6 * h(j, 1)) + _
                   (y(j + 1, 1) / h(j, 1) - h(j, 1) / 6 * z(j + 1, 1)) * (xx(i, 1) - x(j, 1)) + _
                   (y(j, 1) / h(j, 1) - h(j, 1) / 6 * z(j, 1)) * (x(j + 1, 1) - xx(i, 1))
        'i = i + 1
    Next

    NaturalCubicSpline = yy
End Function

'------------------------------------------------------------------
' Procedure : GaussEliminationWithPivoting
' DateTime  : 10.09.2008 13:25
' Author    : Tobias - tobiasschmid.de - vba-blog.de
' Purpose   : Solving an system of equations using GaussElimination and Pivoting
'------------------------------------------------------------------
Public Function GaussEliminationWithPivoting(Matrix() As Variant, Vektor() As VariantAs Variant
    Dim L() As Double    'lower diagonal
    Dim U() As Double    'upper diagonal
    Dim y() As Double
    Dim x() As Double
    Dim Row() As Long, Column() As Long
    Dim k As Long, i As Long, j As Long
    Dim Lij As Double
    Dim temp As Double
    Dim dblMax As Double, rowMax As Long, columnMax As Long

    'only valid for symmetric matrix
    If UBound(Matrix, 1) <> UBound(Matrix, 2) Then Exit Function

    'initialize
    ReDim L(UBound(Matrix, 1), UBound(Matrix, 1))
    ReDim U(UBound(Matrix, 1), UBound(Matrix, 1))
    ReDim Row(UBound(Matrix, 1))
    ReDim Column(UBound(Matrix, 1))

    'GaussElimination method:
    'A  x  = b
    'L U x = b
    'L  y  = b   (Lij = aij/ajj, current Matrix)
    'U  x  = y

    'initialize index-vectors, upper- & lower-diagonal matrix
    For i = 1 To UBound(Matrix, 1)
        'Initialize Index-Vectors
        Row(i) = i
        Column(i) = i
        For j = 1 To UBound(Matrix, 2)
            'initialize matrix U
            U(i, j) = Matrix(i, j)
            'initialize matrix L
            If i = j Then
                L(i, j) = 1
            Else
                L(i, j) = 0
            End If
        Next
    Next

    'get upper- & lower-diagonal matrix
    For k = 1 To UBound(Matrix, 1) - 1    'k: loop over all rows
        'Search for biggest entry in current row and column
        dblMax = Abs(U(Row(k), Column(k)))
        rowMax = k
        columnMax = k
        For i = k + 1 To UBound(Matrix, 1)
            If Abs(U(Row(i), Column(k))) > dblMax Then
                dblMax = Abs(U(Row(i), Column(k)))
                rowMax = i
                columnMax = k
            End If
            If Abs(U(Row(k), Column(i))) > dblMax Then
                dblMax = U(Row(k), Column(i))
                rowMax = k
                columnMax = i
            End If
        Next
        'swap rows and columns
        temp = Row(k)
        Row(k) = rowMax
        Row(rowMax) = temp
        temp = Column(k)
        Column(k) = columnMax
        Column(columnMax) = temp
        'get upper- & lower-diagonal matrix
        For i = k + 1 To UBound(Matrix, 1)   'i: row index
            Lij = U(Row(i), Column(k)) / U(Row(k), Column(k))    'U(i,k) = most left non zero entry in line i
            For j = 1 To UBound(Matrix, 2)   '
                U(Row(i), Column(j)) = U(Row(i), Column(j)) - Lij * U(Row(k), Column(j))
            Next
            L(Row(i), Column(k)) = Lij
        Next
    Next
    'Gauss-Elimination finished

    'get y
    'L  y  = b
    '1 0 0       y   b
    'L 1 0 times y = b
    'L L 1       y   b
    ReDim y(UBound(Vektor, 1), 1)
    For i = 1 To UBound(y, 1)
        temp = 0
        'y(i, 1) = (Vektor(i, 1)
        For j = 1 To i - 1
            temp = temp + y(Row(j), 1) * L(Row(i), Column(j))
        Next
        y(Row(i), 1) = (Vektor(Row(i), 1) - temp)
    Next

    'get x
    'U  x  = y
    'u u u       x   y
    '0 u u times x = y
    '0 0 u       x   y
    ReDim x(UBound(Vektor, 1), 1)
    For i = UBound(x, 1) To 1 Step -1
        temp = 0
        For j = UBound(Vektor, 1) To i + 1 Step -1
            temp = temp + x(Row(j), 1) * U(Row(i), Column(j))
        Next
        x(Row(i), 1) = (y(Row(i), 1) - temp) / U(Row(i), Column(i))
    Next

    'return Result
    For i = 1 To UBound(Vektor, 1)
        'y is temporary vector
        y(i, 1) = x(Row(i), 1)
    Next
    GaussEliminationWithPivoting = y
End Function

Code eingefügt mit Syntaxhighlighter 4.15
AnhangGröße
NaturalCubicSpline.jpg97.91 KB
NaturalCubicSpline.xls56 KB
NaturalCubicSpline.bas7.02 KB
Your rating: Keine Average: 4.3 (6 votes)