Function P2C determines the coefficients of the polynomial which contains N given points. The polynomial is degree N-1. The attached P2C.ods document contains the source for the function and an example of how it is used to find the polynomial Y = 2X⁴ – 6X³ – 8.5X² + 37.5X – 25 through the five points (-0.5,-45) (0,-25) (0.5,-9) (3,11) (3.5,45). The source for the function is also shown below. In the attachment, cells D2:D6 contain the array formula =P2C(A2:A6;B2:B6). Cells A2:A6 have the five X coordinates of the points and B2:B6 have the Y coordinates.
If you want to execute the P2C function in the attachment, you must ensure that option you have set for OpenOffice → Security → Macro Security allows the function in the attachment to run. Press the Help button the the dialog for more information. Options are set with OpenOffice → Preferences on a Mac, Tools → Options on other platforms. You will probably want to copy the the Basic statements to My Macros so you can use the function in any spreadsheet.
Getting started with macros
Rem Purpose: Determine coefficients of polynomial passing through N points
Rem Syntax: P2C(Xvalues,Yvalues)
Rem Xvalues: Range providing X values of the N points
Rem Yvalues: Range providing Y values of the N points
Rem Number of Y values must be equal to number of X values
Rem X values and Y values must be ranges of shape N by 1 or shape 1 by N
Rem All X values must be distinct
Rem Returns: Array of coefficients in ascending order of exponent, constant first
Rem Shape of returned array matches the shape of the Y values
Rem Sample: P2C({-1;1;2};{17;7;11}) is {9;-5;3}
Rem Source: https://en.wikipedia.org/wiki/Divided_differences
Rem Method: For an polynomial of degree i, the coefficient of highest term Xⁱ will be given
Rem by [Y₀,Y₁,…,Yᵢ]. Call this value Qᵢ. Consider the reduced polynomial formed
Rem by dropping the highest term Qᵢ•Xⁱ. Now the highest term is a multiple of Xⁱ⁻¹.
Rem For each of the given points (xᵢ,yᵢ), calculate a new point on the reduced
Rem polynomial as (xᵢ,yᵢ-Qᵢ•xᵢ). One of the points is redundant. This program drops
Rem the last of them. The process is repeated with the reduced polynomial until a
Rem single point (x₀,c) remains representing a polynomial of the form Y=c. This gives
Rem the constant term of the desired polynomial, the other Qᵢ being determined earlier.
Rem Example: The polynomial to be discovered is Y=3X²-5X+9 and we are given points (-1,17)
Rem (1,7) (2,11). x₀=-1 x₁=1 x₂=2, y₀=17 y₁=7 y₂=22. Divided differences give the
Rem highest (X²) coefficient. [y₀,y₁] = (y₁-y₀)/(x₁-x₀) = -5, [y₁,y₂] = (y₂-y₁)/(x₂-x₁)
Rem = 4, [y₀,y₁,y₂] = ([y₁,y₂]-[y₀,y₁])/(x₂-x₀) = 3, the coefficient for X². 3X²
Rem subtracted from (-1,17) and (1,7) gives new points (-1,14) and (1,4), which will
Rem be on a linear polynomial. Divided differences on these points gives the next
Rem (X¹) coefficient. [y₀,y₁] = (y₁-y₀)/(x₁-x₀) = -5, the coefficient for X¹.
Rem Subtracting -5X from (-1,14) gives point (-1,9) so the constant term (X⁰) is 9.
Rem The polynomial coefficients for X⁰,X¹,X² are 9,-5,3.
Option Explicit ' Variables must be declared
Dim NX As Double, NY As Double ' Sizes of X and Y arrays
Dim X0() As Double, Y0() As Double ' 1-dimensional X and Y arrays
Dim DD() As Double ' Divided difference array
Dim I As Integer, J As Integer, K As Integer ' Loop counters
Const E00 = "Don't pass a rectangular range to function P2C" ' Return error message
Sub GetValues(_In() As Variant,Out() As Double,Count As Double) ' Create output array
Select Case 1 ' Determine array dimension
Case UBOUND(_In,1) ' 1 row, N columns
Count = UBOUND(_In,2) ' Size of 1-based array
ReDim Out(Count-1) As Double ' 1-dimension, 0 to Count-1
For K=1 To Count : Out(K-1) = _In(1,K) : Next K ' Copy values, offset 0
Case UBOUND(_In,2) ' N rows, 1 column
Count = UBOUND(_In,1) ' Size of 1-based array
ReDim Out(Count-1) As Double ' 1-dimension, 0 to Count-1
For K=1 To Count : Out(K-1) = _In(K,1) : Next K ' Copy value, offset 0
Case Else ' M rows, N columns
Count = 0 ' Indicate error
End Select ' End Select statement
End Sub ' End GetValues
Function P2C(X1 As Variant,Y1 As Variant) As Variant ' Called with X and Y values
Call GetValues(X1,X0,NX) ' X0 is 0-based array
If NX=0 Then P2C = E00 : Exit Function ' Range was passed if NX is 0
Call GetValues(Y1,Y0,NY) ' Y0 is 0-based array
If NY=0 Then P2C = E00 : Exit Function ' Range was passed if NY is 0
If NX<>NY Then ' Sizes don't agree
P2C = NX & " X values and " & NY & " Y values for function P2C" ' Return error message
Exit Function ' Leave immediately
End If ' End of agreement test
ReDim DD(NX-1,NX-1) As Double ' Divided differences (DD)
For I=NX-1 To 0 Step -1 ' Loop for each polynomial
For K=0 To I : DD(0,K) = Y0(K) : Next K ' DD level zero is Y value
For J=1 To I ' Loop for remaining levels
For K=0 To I-J ' Loop for DD item
DD(J,K) = (DD(J-1,K+1)-DD(J-1,K))/(X0(K+J)-X0(K)) ' Divided difference
Next K ' End item loop
Next J ' End level loop
Y0(I) = DD(I,0) ' Save coefficient
For J=0 To I-1 : Y0(J) = Y0(J)-Y0(I)*X0(J)^I : Next J ' Points on reduced polynomial
Next I ' End polynomial loop
If UBOUND(Y1,1)=1 Then ' Was Y array horizontal?
For K=1 To UBOUND(Y1,2) : Y1(1,K) = Y0(K-1) : Next K ' 1 row, N columns
Else ' No, Y array was vertical
For K=1 To UBOUND(Y1,1) : Y1(K,1) = Y0(K-1) : Next K ' N rows, 1 column
End If ' End direction test
P2C = Y1 ' Return the coefficients
End Function ' End P2C
OpenOffice.org Macros Explained, page 514
Passing arguments to a macro
This topic is intended to illustrate how to pass values to and from a Basic function. For curve fitting with more than a few points, more useful results are typically obtained with Spline interpolation than with polynomial interpolation.