Roger Hui wrote:
> Quite often the best "explanation" for a sequence
> is not a low degree polynomial. An example from
> several days ago was s=: [: p: _1 + 2 ^ i.
> s1=: 2^i. would be a even shorter example.
Taking Roger's advice to heart, we may as well know how to rule out or take
advantage of a "polynomial" sequence, see below.
One kind of polynomial interpolation uses the Newton series, a polynomial
discussed on pages 190-191 of Graham, Knuth, and Patashnik, Concrete
Mathematics
Second Edition. We use a Newton series to investigate the sequence that begins
]data =: 3 ^~ i. 11
0 1 8 27 64 125 216 343 512 729 1000
6 difftb data NB. successive differences of first column
0 1 6 6 0 0
1 7 12 6 0 0
8 19 18 6 0 0
27 37 24 0 0 0
64 61 0 0 0 0
125 0 0 0 0 0
As the third differences are constant, the sequence comes from a third degree
polynomial.
p =: (3 interp data)"0 NB. conjunction interp creates third degree polynomial
data ,: p i. 11 NB. polynomial verb reproduces sequence
0 1 8 27 64 125 216 343 512 729 1000
0 1 8 27 64 125 216 343 512 729 1000
p 0.9 NB. and interpolates between sequence values
0.729
The algebra formula for (p x) can be found using the first row of the
difference
table. In algebra notation (p x) is
p(x) = 0 + 1 x + 6 (x*(x-1))/(1*2) + 6 (x*(x-1)*(x-2))/(1*2*3)
If you expand this and simplify you should get
p(x) = x^3
Here are the verbs and conjunction.
diff =: 2&(-~/\) NB. differences of a vector
difftb =: 4 : '|: diff^:(i.x) x {. y' NB. difference table, x by x
NB. y is data vector with x <: # y .
interp =: 2 : '(y !~ i. >: m) +/ . * {. (>: m) difftb n'
NB. usage m interp data creates a polynomial verb to approximate the
NB. data, a list. If the data is generated by a polynomial of degree
NB. m or less, polynomial p =: m interp data reproduces the data
NB. exactly. Otherwise it approximates the data.
If you look closely at the definition of interp you will see how the formula
for
p(x) shown above is used by interp to create the polynomial
p =: (m interp data)"0. It is best to give the polynomial rank 0 as I just did.
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
