Repost: Cubic Spline+its integral and derivatives
David J. Braden
=======
Over the last few months people have asked about more general interpolation
schemes, as well as "what's the derivative at x of the function F that Excel
uses when smoothing?". It seems Excel uses a cubic spline for interior points,
but dampens it at the edges of the range. In spite of a *lot* of phone calls
and e-mail, no one at Microsoft was able/willing to tell me the algorithm they
use for smoothing data for XY charts. The code below gives you a true cubic
spline of your data, as well as the first three derivatives (the fourth and
beyond are, of course, 0). There is also a function for integrating the cubic spline.
I know it seems like a lot of stuff, but actually it isn't that much: a large
portion is there just for validating input. There's a
compiler directive for shutting off error checking if you wish.
To use:
Simply copy the code into a VBA module. See its comments for an example.
Note: the second derivative of the function is not a differentiable function at
all points, so if you specify points to evaluate that are at "kinky" spots then
you get an #NA for the third derivative.
Error-checking seems to be very tight. If you encounter any bugs or weak spots,
please let me know.
I'm sorry for not having tested this for robustness in international usage;
I don't have the tools to do so. The math seems absolutely sound. If you do
run into a problem in an international setting, please let me know. Except for
the error messages, this will be a non-issue with Excel 2000 and beyond.
Enjoy!
Dave Braden
=========================================
Option Explicit
'This module contains routines for cubic spline interpolation and integration.
'Designed for Microsoft Excel 97 and beyond
'Written 1999/6/30, David J. Braden
'Revisions:
' 1999/7/2 DJB Made 2nd derivative at upper endpoint exact
' Tightened up input validation code
' 1999/7/6 DJB Streamlined (optional) error-checking
'About 20% of the code is for checking that input is valid, hence the length.
'If you "know" that input is always valid, set the compiler directive
' fValidateInput = True.
'To test:
'Into cells A3:A9, enter 1,2,3,4,5,6,7, and into B3:B9, enter 3,3,1,2,1,3,3.
'Into cell B1 enter =6/49 (we will go from 1 to 7 in 50 steps). Into Cell C1,
'enter 1; into C2, enter =C1+$B$1, and drag it down to C50. You should end up
'with 7. While there, go ahead and enter 7 into the last cell (to replace the
'proximal value). Select D1:G50 (that's right, 50 rows, 4 columns), in the
'formula bar type =SplineData(A3:B9,C1:C50), and do a control-shift-enter. The
'first column is a cubic spline interpolation of your data; each subsequant
'column is a higher-order derivative. Of course you don't need to return so much
'data. If you want the cubic-spline interpolation at a single point, no problem:
'specify a single point. If you want more information about the behavior of the
'fitted function evaluated at that point, select up to 4 contiguous cells in the
'same row, and go for it. If you want to know some info about one of the
'derivatives of such a fit, INDEX is very useful.
'For the integral, enter =SplineIntegrate(A3:B9) into a cell. You have the option
'of specifying limits, of course.
'Both functions have optional specification of end-point first derivatives.
#Const fReportErrors = False
#Const fValidateInput = True
Const UseNatBound = 1E+30
Const sTtl = "Spline Evaluation"
Public Function SplineData(XYtbl As Range, EvalPts As Variant, _
Optional EndPtDeriv As Boolean = False, _
Optional LowFirstDeriv As Double = UseNatBound, _
Optional HighFirstDeriv As Double = UseNatBound) As Variant
'This Evaluates a vector of points EvalPts at the cubic spline interpolation of
'the second column of XYtbl to its first column.
'Returns:
' An n by 4 table, where n is the number of EvalPts. Columns contain
' F(x), and the first three derivatives of F wrt x, evaluated at x.
'Optional paramters:
' EndPtDeriv:
' If set to true, will calculate (directional) derivatives at x on the
boundary of the range.
' Default setting is to not do so.
' LowFirstDeriv, HighFirstDeriv:
' Specifications to derive spline such that F'(x) has desired valued at
low/hi end of table.
Const EPS = 0.000000000001 'used to check for kinkiness of f'' when deriving f'''
Dim xy As Variant, ret() As Variant, xEval As Variant, var As Variant
Dim dDel As Double, dA As Double, dB As Double, dY2() As Double
Dim i As Integer, j As Integer, lXYrows As Long
Dim lo As Long, n As Long
Dim fEndpoint As Boolean
xy = XYtbl
If BadXY(xy, XYtbl.Address, lXYrows) Then SplineData = CVErr(xlErrValue): Exit Function
'EvalPts must be a scalar or a vector, to be mapped into an (n,1) table
Vectorize EvalPts, xEval, n
#If fReportErrors Then
If n = 0 Then
MsgBox "The second argument to SplineDataEval" & vbCr & _
"must be a scalar or a vector.", vbExclamation, sTtl
End If
#End If
ReDim dY2(1 To lXYrows) As Double 'these are the 2nd derivatives
ReDim ret(1 To n, 0 To 3) As Variant
spline xy, LowFirstDeriv, HighFirstDeriv, dY2()
For i = 1 To n
var = xEval(i, 1)
If Not Application.IsNumber(var) Then
For j = 0 To 3: ret(i, j) = CVErr(xlErrValue): Next
ElseIf var < xy(1, 1) Or var > xy(lXYrows, 1) Then
For j = 0 To 3: ret(i, j) = CVErr(xlErrNum): Next
Else
lo = Application.Match(var, XYtbl.Columns(1), 1)
fEndpoint = (lo = lXYrows)
If fEndpoint Then lo = lo - 1 'at upper bound of table's domain
dDel = xy(lo + 1, 1) - xy(lo, 1)
dA = (xy(lo + 1, 1) - var) / dDel
dB = 1 - dA
' calculate y(x)
If fEndpoint Then
ret(i, 0) = xy(lXYrows, 2)
Else
ret(i, 0) = dA * xy(lo, 2) + dB * xy(lo + 1, 2) _
+ ((dA * dA * dA - dA) * dY2(lo) _
+ (dB * dB * dB - dB) * dY2(lo + 1)) * dDel * dDel / 6#
End If
fEndpoint = fEndpoint Or (lo = 1 And dB = 0)
If fEndpoint And Not (EndPtDeriv) Then
For j = 1 To 3: ret(i, j) = CVErr(xlErrNA): Next
Else 'have interior point, or will accept a one-sided derivative
'1st derivative
ret(i, 1) = (xy(lo + 1, 2) - xy(lo, 2)) / dDel + dDel / 6# _
* ((3# * dB * dB - 1) * dY2(lo + 1) - (3# * dA * dA - 1#) * dY2(lo))
'2nd derivative
ret(i, 2) = dA * dY2(lo) + dB * dY2(lo + 1)
'3rd derivative
If fEndpoint Then
ret(i, 3) = CVErr(xlErrNA)
Else
dA = (dY2(lo + 1) - dY2(lo)) / dDel 'dA used as temp here
' if in interior, or on a smooth boundary point, all ok
If var > xy(lo, 1) Then
ret(i, 3) = dA
ElseIf Abs(dA - (dY2(lo) - dY2(lo - 1)) / dDel) < EPS Then
ret(i, 3) = dA
Else
ret(i, 3) = CVErr(xlErrNA)
End If
End If
End If
End If
Next
SplineData = ret
End Function
Public Function SplineIntegrate(XYtbl As Range, _
Optional LowX As Double = -1E+306, Optional HighX As Double = 1E+306, _
Optional LowFirstDeriv As Double = UseNatBound, _
Optional HighFirstDeriv As Double = UseNatBound) As Variant
'This integrates the cubic spline interpolation of the second column of XYtbl
'to its first column.
'Returns:
' Integral of fitted function, evaluated from LowX to HighX.
'Optional paramters:
' LowX, HighX:
' Region over which to integrate. If LowX (HighX) isn't specified, then
the lowest (highest)
' point of the domain will be used.
' LowFirstDeriv, HighFirstDeriv:
' Specifications to derive spline such that F' has desired valued at
low/hi end of table.
Dim xy As Variant
Dim dY2() As Double, dA As Double, ret As Double
Dim n As Long, lXYrows As Long, LowXpart As Long, HighXpart As Long
Dim fReverseInt As Boolean
xy = XYtbl
If BadXY(xy, XYtbl.Address, lXYrows) Then
SplineIntegrate = CVErr(xlErrValue): Exit Function
End If
If LowX > HighX Then
dA = HighX: HighX = LowX: LowX = dA: fReverseInt = True
End If
If HighX < xy(1, 1) Or LowX > xy(lXYrows, 1) Then
SplineIntegrate = CVErr(xlErrNum): Exit Function
End If
If LowX < xy(1, 1) Then LowX = xy(1, 1)
If HighX > xy(lXYrows, 1) Then HighX = xy(lXYrows, 1)
ReDim dY2(1 To lXYrows) As Double 'these are the second derivatives
spline xy, LowFirstDeriv, HighFirstDeriv, dY2()
LowXpart = Application.Match(LowX, XYtbl.Columns(1), 1)
HighXpart = Application.Match(HighX, XYtbl.Columns(1), 1)
If HighXpart = lXYrows Then HighXpart = HighXpart - 1
If LowXpart = HighXpart Then 'points are in the same partition
ret = IntEval(HighX, xy, dY2, HighXpart) - IntEval(LowX, xy, dY2, LowXpart)
Else
'integrate from LowX to its upper partition boundary. and
'add analogous qty for HighX
ret = IntEval(HighX, xy, dY2, HighXpart) - IntEval(xy(HighXpart, 1),
xy, dY2, HighXpart) _
+ IntEval(xy(LowXpart + 1, 1), xy, dY2, LowXpart) - IntEval(LowX,
xy, dY2, LowXpart)
For n = LowXpart + 1 To HighXpart - 1
dA = xy(n + 1, 1) - xy(n, 1)
ret = ret + dA * (12 * (xy(n + 1, 2) + xy(n, 2)) - dA * dA * (dY2(n)
+ dY2(n + 1))) / 24#
Next
End If
If fReverseInt Then ret = -ret
SplineIntegrate = ret
End Function
Private Sub Vectorize(vSrc As Variant, x As Variant, lngLen As Long)
Dim intCol As Integer
lngLen = 0
On Error Resume Next
x = vSrc
lngLen = UBound(x, 1): intCol = UBound(x, 2)
If intCol = 0 Then
If lngLen = 0 Then
ReDim x(1 To 1, 1 To 1)
x(1, 1) = vSrc
lngLen = 1
Else
x = Application.Transpose(x)
lngLen = intCol
End If
ElseIf intCol > 1 Then
If lngLen > 1 Then
lngLen = 0
Else
x = Application.Transpose(x)
lngLen = intCol
End If
End If
End Sub
Private Sub spline(xy As Variant, yp1 As Double, ypn As Double, y2() As Double)
' VBA implementation of that found on pp 115-116,
' "Numerical Recipes in C", 2nd ed., by Press et al., Cambridge
University Press
' 1999/6/25 David J. Braden
'xy is column-oriented
Dim n As Long, i As Long, k As Long
Dim p As Double, qn As Double, sig As Double, un As Double, u() As Double
n = UBound(xy)
ReDim u(1 To n - 1) As Double
If yp1 >= UseNatBound Then
y2(1) = 0#: u(1) = 0#
Else
y2(1) = -0.5
u(1) = (3# / (xy(2, 1) - xy(1, 1))) * ((xy(2, 2) - xy(1, 2)) _
/ (xy(2, 1) - xy(1, 1)) - yp1)
End If
For i = 2 To n - 1
sig = (xy(i, 1) - xy(i - 1, 1)) / (xy(i + 1, 1) - xy(i - 1, 1))
p = sig * y2(i - 1) + 2#
y2(i) = (sig - 1#) / p
u(i) = (xy(i + 1, 2) - xy(i, 2)) / (xy(i + 1, 1) - xy(i, 1)) _
- (xy(i, 2) - xy(i - 1, 2)) / (xy(i, 1) - xy(i - 1, 1))
u(i) = (6# * u(i) / (xy(i + 1, 1) - xy(i - 1, 1)) - sig * u(i - 1)) / p
Next
If ypn >= UseNatBound Then
qn = 0#: un = 0#
Else
qn = 0.5
un = (3# / (xy(n, 1) - xy(n - 1, 1))) * (ypn - (xy(n, 2) _
- xy(n - 1, 2)) / (xy(n, 1) - xy(n - 1, 1)))
End If
y2(n) = (un - qn * u(n - 1)) / (qn * y2(n - 1) + 1#)
For k = n - 1 To 1 Step -1: y2(k) = y2(k) * y2(k + 1) + u(k): Next
End Sub
Private Function IntEval(ByVal x As Variant, ByRef xy As Variant, _
ByRef dY2() As Double, ByRef part As Long)
'Evaluate the indefinite integral of a cubic spline fit of second column of
xy to its
'first column, at x, using vector of 2nd derivatives dY2. "part" is the
partition of
'the table that x is in.
Dim dDelta As Double, dTemp As Double, ret As Double
dDelta = xy(part + 1, 1) - xy(part, 1)
dTemp = xy(part + 1, 1) - x: dTemp = (0.5 * dTemp * dTemp)
ret = xy(part, 2) * dTemp / dDelta + dY2(part) _
* (dTemp * dTemp / dDelta - dTemp * dDelta) / 6#
dTemp = x - xy(part, 1): dTemp = (0.5 * dTemp * dTemp)
ret = xy(part + 1, 2) * dTemp / dDelta + dY2(part + 1) _
* (dTemp * dTemp / dDelta - dTemp * dDelta) / 6# - ret
IntEval = ret
End Function
Private Function BadXY(xy As Variant, sXYaddr As String, lngLen As Long) As Boolean
Dim n As Long, intCol As Integer
On Error Resume Next
lngLen = UBound(xy, 1): intCol = UBound(xy, 2)
#If fValidateInput Then
If intCol <> 2 Or lngLen < 2 Then
#If fReportErrors Then
MsgBox "The data in " & sXYaddr & " must be in a contiguous range" &
vbCr & _
"two columns wide, and at least two rows deep." & vbCr & _
"X-values (first column) must all be coercible to numbers.",
vbExclamation, sTtl
#End If
BadXY = True
Else
'make sure x-values in table are positive monotone
For n = 2 To lngLen
If xy(n, 1) <= xy(n - 1, 1) Then Exit For
Next
If n <= lngLen Then
#If fReportErrors Then
MsgBox "The X-values (first column) of " & sXYaddr & _
" must be strictly ascending.", vbExclamation, sTtl
#End If
BadXY = True
End If
End If
#End If
End Function
Chip Pearson
"David J. Braden" <bra...@rochester.rr.com> wrote in message
news:3A11AD33...@rochester.rr.com...
> I repost this by public request. --- DJB
Is this THE David J Braden?? Welcome back to fraternity, Dave. You've been
missed. I hope all is well in your corner of the world. Give me a call
sometime.
--
Cordially,
Chip Pearson
Microsoft MVP - Excel
Pearson Software Consulting, LLC
www.cpearson.com ch...@cpearson.com
Charles Williams
By chance I just found Msoft Excel Appnote Q103841 Generating Smooth curves
in Charts which points to a download (WE0820) and includes references to the
algorithm used. However this is very old and refers to Excel 4.
Maybe this is the same algorithm as they currently use?
regds
Charles
_____________________________________
Decision Models Ltd,17 Binswood Avenue,
Leamington Spa,Warks CV32 5SE, UK
Tel: (44)01926-334289 Fax: (44)01926-881487
David J. Braden <bra...@rochester.rr.com> wrote in message
news:3A11AD33...@rochester.rr.com...