Birkhoff Interpolation Problem

birkoff [1]
Laureano Gonzalez-Vega gvega@matesco.unican.es (Feb. 1998)

The hypersurfaces come from the so called "Birkhoff Interpolation Problem" and in general the problem is to decide if a determinant of a polynomial matrix (in x0,x1,...,xn, for the univariate case) can have, or not, a real point verifiying oftenly some condition (x0<x1<...<xn, for the univariate case).

Here you have one example:

Here EE denotes the information known about the the function (there is one row for every node and one column for every derivative of the polynomial to compute). In this case we look for a degree 7 polynomial
> EE:=array(1..5,1..10,[[1,0,1,0,0,0,0,0,0,0],[0,0,1,0,0,0,0,0,0,0],
                        [1,0,1,0,0,0,0,0,0,0],[0,0,1,1,0,0,0,0,0,0],
                        [1,0,0,0,0,0,0,0,0,0]]);

                  [1    0    1    0    0    0    0    0    0    0]
                  [                                              ]
                  [0    0    1    0    0    0    0    0    0    0]
                  [                                              ]
            EE := [1    0    1    0    0    0    0    0    0    0]
                  [                                              ]
                  [0    0    1    1    0    0    0    0    0    0]
                  [                                              ]
                  [1    0    0    0    0    0    0    0    0    0]
This is the polynomial that if has not real zeros (with ti¹ 0) then the interpolation scheme given by EE is GOOD (always one and only one solution). In this case ti=xi+1-xi.

To apply the code (appended below), it is important that the number of 1's of EE must be equal to the degree + 1.

> FinalPolynomial(EE,[x1,x2,x3,x4,x5],[a0,a1,a2,a3,a4,a5,a6,a7],
                     [t1,t2,t3,t4]);

    4   4   2       2   8       8   2       2   4   2   2       2   2   2   4
2 t1  t3  t4  + 4 t1  t3  + 3 t2  t4  - 4 t2  t1  t3  t4  - 3 t2  t1  t3  t4

            4   2   2   2       4   2   4     2   2   6       4   2   4
     + 10 t2  t1  t3  t4  - 5 t2  t3  t4  - t2  t1  t4  - 2 t2  t1  t4

            2   2   6       4   6       2   6   2        2   4   4
     + 16 t2  t3  t4  + 3 t2  t4  + 2 t2  t3  t4  - 11 t2  t3  t4

           6   2   2       4   4   2       2   4   4        2   8       6   4
     + 8 t2  t3  t4  + 2 t2  t3  t4  + 5 t2  t1  t4  + 18 t2  t4  - 3 t2  t4

            4   6       6   2   2       6   2   2       8   2       6   4
     + 10 t2  t3  - 9 t1  t2  t4  - 6 t1  t3  t4  - 4 t1  t4  + 4 t1  t4

            10       2   6   2       2   6   2       2   4   4       2   6   2
     + 10 t4   + 6 t1  t3  t2  + 2 t1  t3  t4  - 4 t1  t3  t4  + 5 t1  t2  t4

           4   6       10       2   2   6     8   2       6   4       4   6
     + 9 t3  t4  + 3 t3   + 2 t3  t1  t4  + t3  t4  - 5 t3  t4  - 4 t1  t4

           4   2   4        2   8        6   4       4   2   4       2   8
     + 2 t1  t3  t4  + 11 t2  t3  - 10 t2  t3  - 2 t1  t2  t3  + 4 t1  t4

            2   8       4   4   2       2   4   2   2       4   6       10
     + 22 t3  t4  - 3 t1  t2  t4  + 4 t1  t3  t2  t4  - 4 t1  t3  + 4 t1

            8   2        6   4       4   6       2   8       10
     + 13 t1  t2  + 12 t1  t2  - 2 t1  t2  - 8 t1  t2  - 3 t2

            2   6   2        6   2   2       4   4   2       6   4
     - 23 t1  t2  t3  + 19 t1  t2  t3  - 3 t1  t2  t3  + 4 t1  t3

            8   2        2   4   4        8   2
     + 10 t1  t3  - 16 t1  t2  t3  - 11 t2  t3 
The concrete references were to find how to deal with the problem are [2], [3].




Example 1:
with(linalg):
deriv:=proc(P,x,k); if k=0 then P else diff(P,x$k) fi; end:
indexL:=proc(Q,Nodes,i,j)
    local ui,uj;
    ui:=op(i,Nodes);uj:=op(j,Nodes);
    if subs(ui=uj,Q)=0 then 1+indexL(simplify(Q/(ui-uj)),Nodes,i,j)
    else 0; fi; end:
Indicator:= proc(EE,Nodes,Coe)
    local i,j,n,r,x,P,L,M,D;
    n:=nops(Nodes);
    r:=nops(Coe);
    P:= sum('Coe[k]*x**(k-1)', 'k'=1..r);
    L:=[];
    for i from 1 to n do
      for j from 1 to r do
          if EE[i,j]=1 then L:=[op(L),
                subs(x=op(i,Nodes), deriv(P,x,j-1))]  fi; od; od;
    M:=genmatrix(L,Coe);D:=det(M);
    factor(D);
end:
MatrixIndicator:= proc(EE,Nodes,Coe)
    local i,j,n,r,x,P,L,M,D;
    n:=nops(Nodes); r:=nops(Coe);
    P:= sum('Coe[k]*x**(k-1)', 'k'=1..r); L:=[];
    for i from 1 to n do
      for j from 1 to r do
          if EE[i,j]=1 then L:=[op(L),
                subs(x=op(i,Nodes), deriv(P,x,j-1))]  fi; od; od;
    M:=genmatrix(L,Coe);evalm(M);
end:

PoisedIndicator:= proc(EE,Nodes,Coe)
    local lij,ui,uj,i,j,n,r,x,P,L,M,D;
    n:=nops(Nodes); r:=nops(Coe);
    P:= sum('Coe[k]*x**(k-1)', 'k'=1..r); L:=[];
    for i from 1 to n do for j from 1 to r do
          if EE[i,j]=1 then L:=[op(L),
                subs(x=op(i,Nodes), deriv(P,x,j-1))]  fi; od; od;
    M:=genmatrix(L,Coe);D:=det(M);
    if D=0 then RETURN(0); fi;
    for i from 1 to n do
       for j from i+1 to n do
         lij:=indexL(D,Nodes,i,j);
         ui:=op(i,Nodes);uj:=op(j,Nodes);
         D:=simplify(D/(ui-uj)**lij); od; od;
    D/content(D);
end:

FinalPolynomial:= proc(EE,Nodes,Coe,Lengths)
  local i,j,U,U0,T,n;
  U:=PoisedIndicator(EE,Nodes,Coe);
  U0:=subs(op(1,Nodes)=0,U);n:=nops(Nodes);
  for i from 2 to n do
     T:=sum('Lengths[k]**2', 'k'=1..i-1);
     U0:=expand(subs(op(i,Nodes)=T,U0)) od;
  U0; 
end:

# EE:=array(1..5,1..10,[[1,0,1,0,0,0,0,0,0,0],[0,0,1,0,0,0,0,0,0,0],
#                         [1,0,1,0,0,0,0,0,0,0],[0,0,1,1,0,0,0,0,0,0],
#                         [1,0,0,0,0,0,0,0,0,0]]);
# FinalPolynomial(EE,[x1,x2,x3,x4,x5],[a0,a1,a2,a3,a4,a5,a6,a7],
#                     [t1,t2,t3,t4]);

Problem: Decide if the real hypersurface is empty or not.


References

[1]
D. Bini and B. Mourrain. Polynomial test suite. 1996.

[2]
L. Gonzalez-Vega. Applying quantifier elimination to the birkhoff interpolation problem. Journal of Symbolic Computation, 22:83--103, 1996.

[3]
L. Gonzalez-Vega and I. Necula. Quantifier elimination and multivariate birkhoff interpolation. In A. Sydow, editor, 15-th IMACS World Congress on Scientific Computation, Modelling and Applied Mathematics, volume 2, pages 739--744. Wissenschaft and Technik Verlag, 1997.