Re: [Axiom-mail] Another beginner's question: Lagrange interpolation

[Martin Rubey]

Dear Alasdair, Dear Vadim!

Vadim: sorry I haven't sent you the paper by Gemignani yet. It turned out that
I do not have a pdf. If you are still interested, please send me a postal
adress and I send you a collection of papers today.  See below on a tiny
demonstration what FFFG can do.

"Alasdair McAndrew" <address@hidden> writes:

> Just to explore Axiom a little further, I was trying to do a Lagrange
> interpolation on two lists:
> 
> LagrangeInterpolation([1,2,3],[4,5,-6])
> 
> This doesn't work - I get an error message which I (as a beginner) don't
> understand.  So then:
> 
> 1)  How do I do this?

see Bill's mail or below.

> 2) Without asking this mail group, where can I find this information in a
> manner which I can understand?

If you do not want to ask on axiom-mail, you should ask on axiom-math. :-)

Seriously, I think the best and most rapid source of information is 

axiom-math or axiom-mail for "how do I do interpolation, symmetric functions,
                              limits, plot graphs, etc."

axiom-developer for "why doesn't axiom compile on my machine" and "I have a
                     week spare time, and I'm sooo bored, what shall I do".

I'd advise you to spend roughly fifteen minutes 

* browsing HyperDoc, 

* using )what things interpol

* searching the pdf version of the book

* searching the mail archives and the wiki

If you are not successful, send mail.

--  a more scientific answer  -------------------------------------------------

interpolation:

Lagrange interpolation as implemented in Axiom is fine for toy demonstrations.
But if you are interested in interpolation in general, I suggest that you use
my implementation of Bernhard Beckermann and George Labahn's FFFG (Fast
Fraction Free Gaussian) elimination algorithm. It is faster in exact domains
and much more general.  The same file also features (in a separate package)
polynomial interpolation using Newton, albeit in an extremely specialised
setting:

      ++ \spad{newton}(l) returns the interpolating polynomial for the values
      ++ l, where the x-coordinates are assumed to be [1,2,3,...,n] and the
      ++ coefficients of the interpolating polynomial are known to be in the
      ++ domain F. I.e., it is a very streamlined version for a special case of
      ++ interpolation. 

You will have to install my package manually though, or use wh-sandbox. Here
are the docstrings:

    interpolate: (List Fraction D, List Fraction D, NonNegativeInteger) 
                -> Fraction SUP D
      ++ \spad{interpolate(xlist, ylist, deg} returns the rational function with
      ++ numerator degree \spad{deg} that interpolates the given points using
      ++ fraction free arithmetic.

    generalInterpolation: (List D, CoeffAction, 
                           Vector V, List NonNegativeInteger) -> Matrix SUP D
      ++ \spad{generalInterpolation(C, CA, f, eta)} performs Hermite-Pade
      ++ approximation using the given action CA of polynomials on the elements
      ++ of f. The result is guaranteed to be correct up to order
      ++ |eta|-1. Given that eta is a "normal" point, the degrees on the
      ++ diagonal are given by eta. The degrees of column i are in this case
      ++ eta + e.i - [1,1,...,1], where the degree of zero is -1.
      ++
      ++ The first argument C is the list of coefficients c_{k,k} in the 
      ++ expansion <x^k> z g(x) = sum_{i=0}^k c_{k,i} <x^i> g(x).
      ++
      ++ The second argument, CA(k, l, f), should return the coefficient of x^k
      ++ in z^l f(x).


So, polynomial interpolation is interpolate([1,2,3],[4,5,-6], 3), Hermite-Pade
is, for example

(1) -> f :=reduce(+, [monomial(random 40, i)$SUP INT for i in 0..5])

           5      4      3      2
   (1)  12?  + 25?  + 29?  + 13?  + 20? + 1
                                     Type: SparseUnivariatePolynomial Integer
(2) -> F ==> FFFG(INT, SUP INT)
                                                                   Type: Void
(3) -> eta := [2,1,1]

   (3)  [2,1,1]
                                                   Type: List PositiveInteger
(4) -> n := reduce(+, eta) 

   (4)  4
                                                        Type: PositiveInteger
(5) -> M := generalInterpolation(DiffC(n)$F, DiffAction$F, [f, D f, D D f], 
eta)$F

   (5)
   +          2                                                               +
   |- 1045836?  + 939076? + 1012588  3237618? + 3755082   16407212? + 9862196 |
   |                                                                          |
   |            133656               - 1045836? + 299430        1116948       |
   |                                                                          |
   +           - 141758                   - 374757        - 1045836? - 1238506+
                              Type: Matrix SparseUnivariatePolynomial Integer
(6) -> [f, D f, D D f] * M

   (6)
               7            6           5            4
   [- 12550032?  - 14876988?  + 5298712?  + 46971396? ,
             6            5             4
    38851416?  + 63251274?  + 101150172? ,
              6             5             4
    196886544?  + 528526652?  + 538380288? ]
                              Type: Vector SparseUnivariatePolynomial Integer


The last line demonstrates that we have indeed found a solution to the given
(Hermite-Padé) interpolation problem. For example, let [c1, c2, c3] be the
first column of M, then

ord(c1*f + c2 * D f + c3 * D D f) >= 4.


Note that, although seemingly different, interpolate and generalInterpolation
use exactly the same algorithm, with slightly different parameters.


Possible future project:  I have implemented several interpolation algorithms,
each useful in it's setting.  It would be good to design a general user
interface, that is more natural than I currently have.  Unfortunately, I do not
have the time...

Second project: implement Peter Olver's ideas about multivariate interpolation:

Olver, P.J., On multivariate interpolation, Stud. Appl. Math. 116 (2006)
201-240.
http://www.math.umn.edu/~olver/n_/mv.pdf

I'd really like to do this, time permitting.



Martin