Dienstag, 13. Januar 2009

Spline Interpolation

Für die Freunde von Kurven habe ich nachstehend den Quellcode für die Spline Interpolation hinzugefügt. Die Funktion kann ebenfalls, wie die Lineare Interpolation in Excel "eingebaut" werden und wird für Interpolation von Messwerten aus vorhandenen Messwertpaaren benutzt.

Das Programm stammt aus dem Quellcode des Freeware Programms XNUMBERS 5.6 - Multi Precision Floating Point Computing and Numerical Methods for EXCEL.
Ein Muss für jeden Excel-Anwender. Einfach suchen unter: Digilander



====================== ANFANG ======================
  • 'Interpolates one point from sorted X,Y data pairs using Cubic Spline Interpolation
  • Function Interpol_cspline(Xin As Range, Yin As Range, Xtarget As Double)
  • Dim n%, ny% 'Xin and Yin counts
  • Dim i% 'loop counter
  • Dim ihi%, m%
  • Dim a#, b#, p#, qn#, sig#, un#, h0#, h1#, h2#, Yout#
  • ny = Yin.Cells.Count
  • n = Xin.Cells.Count
  • ' Next check to be sure that there are the came counts of Xin and Yin
  • If n <> ny Then
  • Interpol_cspline = CVErr(xlErrRef) 'Uneven counts of Xin and Yin"
  • GoTo err_ret
  • End If

  • ReDim x(1 To n) As Double
  • ReDim Y(1 To n) As Double
  • 'ReDim u(1 To n - 1) As Single
  • ReDim U(1 To n - 1) As Double 'mod. LV 31-8-02
  • ReDim ypp(1 To n) As Double 'these are the 2nd derivative values
  • 'populate the input arrays
  • For i = 1 To n
  • x(i) = Xin(i)
  • If i > 1 Then
  • If x(i) <= x(i - 1) Then 'Check if X values are not repeating and are sorted in ascending order
  • Interpol_cspline = CVErr(xlErrValue) 'Error: Not sorted or repeating
  • GoTo err_ret
  • End If
  • End If
  • Y(i) = Yin(i)
  • Next i
  • ypp(1) = 0 'First knot boundary condition
  • U(1) = 0 'First knot boundary condition
  • For i = 2 To n - 1
  • h0 = x(i + 0) - x(i - 1)
  • h1 = x(i + 1) - x(i - 0)
  • h2 = x(i + 1) - x(i - 1)

  • sig = h0 / h2
  • p = sig * ypp(i - 1) + 2#
  • ypp(i) = (sig - 1) / p
  • U(i) = (Y(i + 1) - Y(i)) / h1 - (Y(i) - Y(i - 1)) / h0
  • U(i) = (6# * U(i) / h2 - sig * U(i - 1)) / p

  • Next i

  • qn = 0
  • un = 0
  • ypp(n) = (un - qn * U(n - 1)) / (qn * ypp(n - 1) + 1) 'Last knot boundary condition
  • For i = n - 1 To 1 Step -1
  • ypp(i) = ypp(i) * ypp(i + 1) + U(i) 'Backfill the 2nd derivatives
  • Next i
  • ''''''''''''''''''''''''''''''''''''''''
  • 'Spline evaluation at target X point
  • '''''''''''''''''''''''''''''''''''''''''
  • 'Find correct interval using halving binary search
  • i = 1
  • ihi = n
  • Do While (i <>
  • m = (i + ihi) \ 2 'Calculate the midpoint
  • If (Xtarget < ihi =" m" i =" m">
  • Loop
  • ' i = beginning of the correct interval
  • h0 = x(i + 1) - x(i) 'Calc the width of the X interval
  • a = (x(i + 1) - Xtarget) / h0
  • b = (Xtarget - x(i)) / h0
  • Yout = a * Y(i) + b * Y(i + 1) + ((a ^ 3 - a) * ypp(i) + (b ^ 3 - b) * ypp(i + 1)) * (h0 ^ 2) / 6#
  • Interpol_cspline = Yout
  • err_ret:
  • End Function
====================== ENDE ======================

Gruß Chiochip

Keine Kommentare: