Previous Page Next Page Contents

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

Option: F

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

Changes




Do you have questions or comments?


Copyright © SciFace Software GmbH & Co. KG 2000