Polynomial Interpolation
Vector 12:2, page 26
Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 26
Polynomial Interpolation¶
Interpolation means the process of using a table of values of an independent variable \(x\) with corresponding function values \(f(x)\), to determine a value of \(f(x)\) at a non-tabular value of \(x\). For many functions, a polynomial which passes through the tabulated points gives an adequate representation near the non-tabular value. Such a polynomial is then used to interpolate values of \(f(x)\) at other values of \(x\) without explicit calculation of its coefficient values.
There is a variety of techniques for doing this, among them Lagrangian polynomials and Neville’s algorithm. The technique used here is known as Newton’s Divided Differences, and consists of taking successive differences of f-values which are divided at each stage by the matching differences in the column of x-values. The values of the nth column of f-differences are divided by the differences of the corresponding x-values which are n apart, for example:
In the 1st diff column 18-6 is divided by 5-1, and 6-18 by 8-5. In the 2nd diff column, 3-2 is divided by 5-0 and -4-3 by 8-1, and so on.
The result of the verb divdifs is x laminated with the values at the head of the columns, viz.
from which the original values of f could, if required, be reconstructed by reversing the arithmetic.
The J realisation which follows is in the style of a Pascal program, which is a natural way of expressing an iterative algorithm of this sort. Indeed this algorithm illustrates the ease with which transition is possible between J and more familiar programming languages.
NB. Obtain leading values in divided difference columns
divdifs =. 0 : 0
:
i =.0 [ f=. a=. y. [ n=.#x. NB. initialization
loop) $.=.>(n > i=.i+1) { end;$. NB. for i = 1 to n-1 do
j =. i-1 NB. start of inner loop
loop1) $. =. >(n > j=.j+1) { end1;$. NB. for j = i to n-1 do
t =. (j{f - (j-1){f)%(j{x. - (j-i){x.) NB. a[j]:=
a =. t j } a NB. ((f[j]-f[j-1])/(x[j]-x[j-i])
$.=.loop1 NB. end of inner loop
end1) f =. a NB. update f
$.=.loop NB. end of outer loop
end) x.,:a NB. result is x,:leading diff'nces.
)
Interpolation is performed by an algorithm which uses the leading divided differences. A merit of this technique is that as the algorithm progresses through its iterations the result vector grows by catenating interpolated values of successively higher polynomial degree as more points are brought into consideration from the right of the data.
NB. Program to perform polynomial interpolations of increasing order
interpol =. 0 : 0
''
:
f =. 1{x. [ x =. 0{x. NB. f = function values, x = x-values
a =. 1{t =. x divdifs f NB. a = leading divided differences
n =. #a [ i =. 0 [ p =. 0{a NB. p =a[0], first estimate
sol =. a NB. sol is result vec (see last NB.)
loop) $. =. >(n>i =.i+1) { end;$. NB. for i = 1 to n-1 do
j =. i [ t =. 1 NB. start of inner loop
loop1) $. =. >(_1<j=.j-1) { end1;$. NB. for j = i-1 downto 0 do
t =. t*(y. -j{x) NB. t := t*(y. -x[j])
$. =. loop1 NB. end of inner loop
end1) p =. p+t*((i){a) NB. p:=p + t*a[i]
sol =. p i } sol NB. sol[i]:=p
$.=.loop NB. end of outer loop
end) sol NB. sol is successive polyn'l interpolns.
)
x =.8 9 9.5 11
f =. 2.079442 2.197225 2.251292 2.397895
xf =. x,:f
xf interpol 9.2
2.079442 2.22078 2.21924 2.21921
The last line shows that the cubic interpolation is 2.21921, although this is little different from the value obtained by quadratic interpolation.
It is not necessary that the x values be in ascending order. For example, using the first example given above the successive polynomial interpolations are
If the points are introduced in x-value order 0 8 1 5, the polynomial interpolations are
which reflects the fact that the straight line approximation will vary greatly according to when the point (5,18) is introduced.
Choleski Decomposition ("Square Root") of a Matrix¶
This is a method of factorisation of a symmetric matrix
where L is a lower triangular matrix. If A is positive definite the elements of L will be real, otherwise some will be complex. The mathematical equations which govern the method are:
These lead to the definition of some subsidiary verbs as follows:
ssqbut1=. +/@( (-:i.@#@]) # *:@] ) NB. 2 ssqbut1 4 2 7 3 5 yields the
NB. sum of squares of 4 2 3 and 5
vb=. %(%:@{.) NB. vb 2 3 4 yields
NB. sqrt(2), 3/sqrt(2), 4/sqrt(2)
stap1=. (vb@{.)&.|: NB. performs vb on first column
NB. of a matrix
column=. (*i.)@# NB. column yields 0 n 2n 3n ...(n-1)n
NB. where n=#y. (y. is a matrix)
fillmatrix =. $~$ NB. x. fillmatrix y. yields matrix of
NB. same size as y. filled with x.
fillM1 =. '' : ' (stap1 y.) (column y.) } (0 fillmatrix y.)'
NB. 1 fillM1 m yields zero
NB. matrix with first column
NB. containing the result of stap1
As with interpol a Pascal program is a natural way of expressing the essentially iterative algorithm, whose development, given the above verbs is now straightforward:
NB. PROGRAM choleski
ch =. 0 : 0
k=. 0 [ n=.#M=. 1 fillM1 y. NB. initialisation
loop1) $. =. >(n>k=.k+1){ end1;$. NB. if k >= n go to end1
t=. %: k{k{y.- (k ssqbut1 (k{M)) NB. t = sqrt(Akk - sum(Mks^2)
M=. t((n+1)*k) } M NB. Mkk:= t
j =. k NB. start of inner loop
loop2) $. =. >(n>j=.j+1) {end2;$. NB. if j>= n go to end2
t=. (j{k{y.- (k{M) +/ .* j{M )% t NB. t:=(Ajk - sum MjkMks)/Mkk
M=. t ((j * n) + k) } M NB. Mjk:=t
$. =. loop2 NB. end of inner loop
end2) NB. j=n
$. =. loop1 NB. end of outer loop
end1)M
:
''
)
Two examples are:
A ]x=.ch A x+/ .* |:x
4 2 14 2 0 0 4 2 14
2 17 _5 1 4 0 2 17 _5
14 _5 83 7 _3 5 14 _5 83
a ch a
1 4 6 1 0 0
4 0 3 4 0j4 0
6 3 2 6 0j5.25 0j2.53722
The authors of this article are post-graduate mathematics students at the Universities of Kaiserslautern, Eindhoven and Strathclyde respectively, and were introduced to J version 7, and to each other, in the course of ECMI (European Consortium for Mathematics in Industry) exchanges at the University of Strathclyde.