Re: Interpolation, bugfix for rationalinterpolation
Subject: Re: Interpolation, bugfix for rationalinterpolation
From: Martin RUBEY
Date: Mon, 7 Oct 2002 12:01:13 +0200 (CEST)
Dear Francisco,
thank you very much for your report, it was a VERY stupid bug, below you
find the corrected version. Note that I changed the syntax of
rationalinterpolation a little, it should be easier to use now... (it is
still not efficient, but I do not have the time right now to implement the
Burlisch-Stoer algorithm) Anyway, now you would call
rationalinterpolation(xvalues, yvalues, x, m, k);
with
xvalues:[-1,0,2];
yvalues:[4,1,-1];
m:2;
k:0;
to get the right result.
PLEASE do change your file, because the results were really only correct
in the case xvalues=[1,2,3,...]
(which was however the only case used by guess...)
/* guess.usg
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* *
* Copyright (C) 2002 Martin Rubey <Martin.Rubey@LaBRI.fr> *
* *
* This file is part of GNU Maxima. *
* *
* This program is free software; you can redistribute it and/or *
* modify it under the terms of the GNU General Public License as *
* published by the Free Software Foundation; either version 2 of *
* the License, or (at your option) any later version. *
* *
* This program is distributed in the hope that it will be *
* useful, but WITHOUT ANY WARRANTY; without even the implied *
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR *
* PURPOSE. See the GNU General Public License for more details. *
* *
* History: *
* This is a translation of the Mathematica package Rate.m *
* by Christian Krattenthaler <Kratt@pap.univie.ac.at>. *
* The translation to Maple was done by Jean-Francois Beraud *
* <Jean-Francois.Beraud@sic.sp2mi.univ-poitiers.fr> and Bruno Gauthier *
* <Bruno.Gauthier@univ-mlv.fr> *
* *
* All features of this package are due to C. Krattenthaler *
* The help text is due to Jean-Francois Beraud and Bruno Gauthier *
* *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
A package to guess closed form for a sequence of numbers.
CALLING SEQUENCE:
guess(l, optional_args);
SYNOPSIS:
- This package provides functions to find a closed form for a
sequence,
of numbers within the hierarchy of expressions of the
form,
<rational function>, <product of rational function>, <product of
product,
of rational function>, etc.
EXAMPLES:
guess([1,2,3]);
[i0]
guess([1,4,9,16]);
2
[i0 ]
guess([1,2,6,24,120]);
i0 - 1
/===\
! !
[ ! ! (i1 + 1)]
! !
i1 = 1
guess(makelist(product(product(GAMMA(i)*i^2,i,1,j),j,1,k),k,1,8));
i0 - 1 i1 - 1 i2 - 1
/===\ /===\ /===\ 2
! ! ! ! ! ! (i3 + 3)
[ ! ! 4 ! ! 18 ! ! ---------]
! ! ! ! ! ! i3 + 2
i1 = 1 i2 = 1 i3 = 1
guess([1,2,7,42,429,7436,218348,10850216]);
i0 - 1 i1 - 1
/===\ /===\
! ! ! ! 3 (3 i2 + 2) (3 i2 + 4)
[ ! ! 2 ! ! -----------------------]
! ! ! ! 4 (2 i2 + 1) (2 i2 + 3)
i1 = 1 i2 = 1
guess(makelist(k^3+k^2,k,1,7));
Dependent equations eliminated: (6)
i0 - 1
/===\
2 ! ! 5040
[i0 (i0 + 1), 2 ! ! (- --------------------------------------),
! ! 4 3 2
i1 = 1 i1 - 24 i1 + 245 i1 - 1422 i1 + 360
i0 - 1
/===\
! ! (i1 + 1) (i1
+ 2)
2 ! !
-----------------]
! ! 2
i1 = 1 i1
Note that the last example produces three solutions. The first and the
last are
equivalent, but the second is different! In this case,
guess(makelist(k^3+k^2,k,1,7),1);
or
guess(makelist(k^3+k^2,k,1,7),"one");
2
find only the solution i0 (i0 + 1), which is a rational function, and is
also
the first function guess finds.
PARAMETERS:
l - a list of numbers,
level - an integer (optional),
"one" - the string "one" (optional),
"nogamma" - the string "nogamma" (optional),
SYNOPSIS:,
- guess(l) tries to find a closed form for a sequence within the
hierarchy,
of expressions of the form <rational function>, <product of
rational,
function>, <product of product of rational function>, etc.
- guess(l,level) does the same thing as guess(l) but it searches only
within the first 'level' levels
- guess(l,"one") does the same thing as guess(l) but it returns the first
solution it finds.
- guess(l,"nogamma") does the same thing as guess(l) but it returns
expressions without GAMMA functions. In fact, there is not much
difference
just at the moment, because Maxima doesn't simplify products yet...
*/
/* guess.mac -*- mode: Maxima; -*-
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* *
* Copyright (C) 2002 Martin Rubey <Martin.Rubey@LaBRI.fr> *
* *
* This file is part of GNU Maxima. *
* *
* This program is free software; you can redistribute it and/or *
* modify it under the terms of the GNU General Public License as *
* published by the Free Software Foundation; either version 2 of *
* the License, or (at your option) any later version. *
* *
* This program is distributed in the hope that it will be *
* useful, but WITHOUT ANY WARRANTY; without even the implied *
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR *
* PURPOSE. See the GNU General Public License for more details. *
* *
* History: *
* This is a translation of the Mathematica package Rate.m *
* by Christian Krattenthaler <Kratt@pap.univie.ac.at>. *
* The translation to Maple was done by Jean-Francois Beraud *
* <Jean-Francois.Beraud@sic.sp2mi.univ-poitiers.fr> and Bruno Gauthier *
* <Bruno.Gauthier@univ-mlv.fr> *
* *
* All features of this package are due to C. Krattenthaler *
* The help text is due to Jean-Francois Beraud and Bruno Gauthier *
* *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
*/
/*
* Rational Interpolation
* Gives the rational interpolant to the set of points defined by xlist
and
* ylist, where m and k are respectively the degrees of the numerator and
* denominator, and xlist is a list of m+k+1 abscissas of the
interpolation
* points, x is the variable the result is supposed to be a function of.
*/
rationalinterpolation(xlist, ylist, x, m, k) :=
block([tempvec : makelist(1, i, 1, m+k+1), /* contains the new row of mat
*/
rowlist, /* first set of rows of mat
*/
mat, /* matrix that describes the interpolation problem
*/
varlist : makelist('x[i], i, 1, m+k+2)],
mode_declare([tempvec,rowlist,varlist,mat],list,[m,k],fixnum),
if max(m, k) > 0
then rowlist : cons(tempvec,
makelist(tempvec : tempvec * xlist, i, 1, max(m,
k)))
else rowlist : [tempvec],
mat : transpose(apply(matrix,
append(rest(rowlist, -(max(m, k) - m) ),
-1 * makelist(rowlist[i] * ylist,
i, 1, k + 1)))),
mat : ev(mat . varlist, SCALARMATRIXP : false),
/* not sure whether it is save to modify xlist... */
xlist : linsolve(makelist(mat[i, 1], i, 1, (m+k)+1), varlist),
if length(xlist) = 0
/* something went wrong */
then NULL
/* use the solution to define the interpolating rational function */
else factor(subst(xlist, sum('x[i+1]*x^i, i, 0, m)
/sum('x[(i+m)+2]*x^i, i, 0, k))));
/* Intermediate function */
guesscons(l, t) :=
block([lsize : length(l), res : [], x, ri],
mode_declare(lsize, fixnum, res, list, ri, any),
for i : 0 thru lsize-2 do
(ri : rationalinterpolation(makelist(k, k, 1, lsize-1), rest(l,-1),
x, (lsize-2)-i, i),
if ri # NULL
then if (subst(x=lsize, denom(ri)) # 0)
and
(subst(x=lsize, ri)-last(l) = 0)
and
not member(subst(x=t, ri), res)
then res : cons(subst(x=t, ri), res)),
res);
/*
* Main function of the package
* it tries to find a closed form for a sequence
* within the hierarchy of expressions of the
* form <rational function>, <product of rational functions>,
* <product of product of rational functions>, etc. It may
* give several answers
*/
guess(l, [optargs]) :=
block([lsize : length(l),
tempres, maxlevel,
maxlevellist : sublist(optargs, numberp),
res : [],
onep : member("one", optargs),
unevp : member("nogamma", optargs), g],
mode_declare([lsize, maxlevel], fixnum,
[tempres, maxlevellist, res], list,
[onep, unevp], boolean, g, any),
optargs : delete("nogamma", delete("one", optargs, 1), 1),
if length(maxlevellist) > 1 or length(optargs)-length(maxlevellist)
> 0
then error("Wrong number of optional arguments: ", optargs)
else maxlevel : mode_identity(fixnum, apply(min, cons(lsize-1,
maxlevellist)) - 1),
g : make_array('ANY, maxlevel + 1),
for k : 0 thru maxlevel do
(g[k] : l,
l : makelist(l[i+1]/l[i], i, 1, (lsize-k)-1),
tempres : guesscons(g[k], concat('i, k)),
if tempres # []
then (if k > 0
then for i : 1 thru k do
if unevp
then tempres :
subst('j = concat('i, (k-i)+1),
map(lambda([expr],
g[k-i][1] *
'product(expr, j, 1,
concat('i,
k-i)-1)),
tempres))
else tempres :
subst('j = concat('i, (k-i)+1),
map(lambda([expr],
g[k-i][1] *
product(expr, j, 1,
concat('i,
k-i)-1)),
tempres)),
res : append(res, tempres),
if onep then return())),
res);