numeric::lagrange
-- polynomial
interpolationnumeric::lagrange
computes an interpolating polynomial
through data over a rectangular grid.
numeric::lagrange(nodes, values, ind <, F>)
nodes |
- | a list [L1,L2,..] of d lists
L.i defining a d-dimensional rectangular grid
{[x1,x2,..]; x1 in L1, x2 in L2, ..}.The lists L.i may have different lengths n.i=nops(L.i). The elements of each L.i must be distinct. |
values |
- | a d-dimensional
array(1..n1,1..n2,..,[...]) associating a value with each
grid point:
[L1[i1],L2[i2],..] --> values[i1,i2,..], i1=1,..,n1, i2=1,..,n2, .. . |
ind |
- | a list of d indeterminates or arithmetical
expressions. Indeterminates are either identifiers (of domain type
DOM_IDENT ) or indexed
identifiers (of type "_index" ). |
|
- | either Expr or any field of
category Cat::Field |
An interpolating polynomial P of domain type DOM_POLY
in the indeterminates specified
by ind
over the coefficient field F
is
returned. The elements in ind
that are not indeterminates
but arithmetical expressions are not used as indeterminates in
P, but enter its coefficients: the polynomial is
``evaluated'' at these points. If no element of ind
is an
indeterminate, then the value of the polynomial at the point specified
by ind
is returned. This is an element of the field
F
or an expression, if F=
Expr.
genpoly
, numeric::cubicSpline
,
poly
ind
=[X1,X2,..]
are specified. The interpolating polynomial P=poly(..,
[X1,X2,..],F) satisfies
evalp(P,X1=L1[i1],X2=L2[i2],..) = value[i1,i2,..]for all points [L1[i1],L2[i2],..] in the grid. P is the polynomial of minimal degree satisfying the interpolation conditions, i.e.,
degree
(P,X.i)<n.i.ind
=[X1,X2,..] and then evaluate
P(v1,v2,..). It is faster to compute this value directly by
numeric::lagrange
with
ind
=[v1,v2,..]. Cf. examples 1 and 3.poly
(.., F)
.F
<>Expr the
grid nodes as well as the entries of values
must be
elements of F
or must be convertible to such elements.
Conversion of the input data to elements of F
is done
automatically.We consider a 1-dimensional interpolation problem. To each node x[i] a value y[i] is associated. The interpolation polynomial P with P(x[i])=y[i] is:
>> L1 := [1, 2, 3]: values := array(1..3, [y1, y2, y3]): P := numeric::lagrange([L1], values, [X])
/ / y1 y3 \ 2 / 5 y1 3 y3 \ poly| | -- - y2 + -- | X + | - ---- + 4 y2 - ---- | X + \ \ 2 2 / \ 2 2 / \ (3 y1 - 3 y2 + y3), [X] | /
The evaluation of P at the point X=5/2 is given by:
>> evalp(P, X = 5/2)
3 y2 y1 3 y3 ---- - -- + ---- 4 8 8
It can also be computed directly without the symbolic polynomial:
>> numeric::lagrange([L1], values, [5/2])
3 y2 y1 3 y3 ---- - -- + ---- 4 8 8
>> delete L1, values, P:
We demonstrate multi-dimensional interpolation. Consider data over the following 2 x 3 grid:
>> XList := [1, 2]: YList := [1, 2, 3]: values := array(1..2, 1..3, [[1, 2, 3], [3, 2, 1]]): P := numeric::lagrange([XList, YList], values, [X, Y])
poly(- 2 X Y + 4 X + 3 Y - 4, [X, Y])
Next, interpolation over a 2 x 3 x 2 grid is demonstrated:
>> L1 := [1, 2]: L2 := [1, 2, 3]: L3 := [1, 2]: values := array(1..2, 1..3, 1..2, [[[1, 4], [1, 2], [3, 3]], [[1, 4], [1, 3], [4, 0]]]): numeric::lagrange([L1, L2, L3], values, [X, Y, Z])
2 2 poly(- 3 X Y Z + 7/2 X Y + 10 X Y Z - 23/2 X Y - 7 X Z + 2 2 8 X + 7/2 Y Z - 3 Y - 27/2 Y Z + 12 Y + 13 Z - 11, [X, Y, Z])
>> delete XList, P, L1, L2, L2, values:
We interpolate data over a 2-dimensional grid:
>> n1 := 4: L1 := [i $ i = 1..n1]: n2 := 5: L2 := [i $ i = 1..n2]: f := (X, Y) -> 1/(1 + X^2 + Y^2): values := array(1..n1, 1..n2, [[f(L1[i], L2[j]) $ j=1..n2] $ i=1..n1]):
First we compute the symbolic polynomial:
>> P := numeric::lagrange([L1, L2], values, [X, Y])
3 4 3 3 poly(- 5563/23108085 X Y + 16376/4621617 X Y - ... - 4401895/3081078 Y + 4199983/2567565, [X, Y])
Fixing the value Y=2.5 this yields a polynomial in X. It can also be computed directly by using an evaluation point for the indeterminate Y:
>> numeric::lagrange([L1, L2], values, [X, 2.5])
3 2 poly(0.0007372500794 X - 0.002155538175 X - 0.03076935248 X + 0.1533997618, [X])
If all indeterminates are replaced by evaluation points, then the corresponding interpolation value is returned:
>> numeric::lagrange([L1, L2], values, [1.2, 2.5])
0.1146465319
>> delete n1, n2, f, values, P:
We demonstrate interpolation over a non-standard field. Consider the following data over a 2 x 3 grid:
>> XList := [3, 4]: YList := [1, 2, 3]: values := array(1..2, 1..3, [[0, 1, 2], [3, 2, 1]]):
With the following call these data are converted to integers modulo 7. Arithmetic over this field is used:
>> F := Dom::IntegerMod(7): P := numeric::lagrange([XList, YList], values, [X, Y], F)
poly(5 X Y + 5 X + 5, [X, Y], Dom::IntegerMod(7))
Evaluation of P at grid points reproduces the associated values converted to the field:
>> evalp(P, X = XList[2], Y = YList[3]) = F(values[2, 3])
1 mod 7 = 1 mod 7
>> delete XList, YList, values, F, P:
{[x1,x2,..]; x1 in L1, x2 in L2, .. }specified by the lists
L.j=[x[j,1],x[j,2],..] of length n.jwith associated values
P(x1[i1], x2[i2], ..) = v[i1,i2,..]the interpolating polynomial in the indeterminates X1,X2,.. is given by
P(X1,X2,..) = sum(sum(.. (v[i1,i2,..]*p[1,i1](X1)*p[2,i2](X2)*..)..))with the Lagrange polynomials
p[j,k](X) = product((X-x[j,l])/(x[j,k]-x[j,l]), l=1..n.j, l<>k)associated with the k-th node of the j-th coordinate.
lagrange