## `numeric::lagrange` -- polynomial interpolation

#### Introduction

`numeric::lagrange` computes an interpolating polynomial through data over a rectangular grid.

#### Call(s)

`numeric::lagrange(nodes, values, ind <, F>)`

#### Parameters

`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"`).

#### Options

``F`` | - | either Expr or any field of category `Cat::Field`

#### Returns

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.

#### Related Functions

`genpoly`, `numeric::cubicSpline`, `poly`

#### Details

- Assume that indeterminates `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.
- If only interpolating values at concrete numerical points X1=v1,X2=v2,.. are required, then we recommend not to compute P with symbolic indeterminates `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.

#### Option: F

- The returned polynomial is of type `poly``(.., F)`.
- For the default field Expr all input data may be arbitrary MuPAD expressions. Standard arithmetic over such expressions is used to compute the polynomial.
- For `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.

#### Example 1

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:

#### Example 2

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:

#### Example 3

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:

#### Example 4

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:

#### Background

- For a d-dimensional rectangular grid

               {[x1,x2,..]; x1 in L1, x2 in L2, .. }

  specified by the lists

      L.j=[x[j,1],x[j,2],..] of length n.j

  with 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.

#### Changes

- used to be called `lagrange`
- Also 1-dimensional grids must now be specified as a list with a list of nodes. The values must now be specified by an array. Also for 1-dimensional interpolation the indeterminate must now be specified by a list.
- Interpolation is now possible over grids of arbitrary dimension. Further, it is now possible to specify indeterminates as well as evaluation points.

Do you have questions or comments?
Copyright © SciFace Software GmbH & Co. KG 2000