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.
- P2C chart.gif (40.06 KiB) Viewed 3573 times
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
- P2C function location.png (107.62 KiB) Viewed 3573 times
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.