(Bugfix) ratinterpol does not work with polynomial expressions
Subject: (Bugfix) ratinterpol does not work with polynomial expressions
From: andre maute
Date: Sat, 28 Dec 2013 06:33:20 +0100
I have added a new interpol.mac, which I believe fixes this issue.
The new version uses setdifference, perhaps there is a better solution.
Andre
------ output with new interpol.mac ------
[user at localhost ~]$ maxima -b ratinterpol-20132712.mac
Maxima 5.31.3 http://maxima.sourceforge.net
using Lisp SBCL 1.1.8-2.fc19
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
STYLE-WARNING: redefining MAXIMA::$FILE_TYPE in DEFUN
(%i1) batch("ratinterpol-20132712.mac")
read and interpret file: /home/user/ratinterpol-20132712.mac
(%i2) display2d:false
(%o2) false
(%i3) load("interpol")
(%o3) "/home/user/opt/maxima/share/maxima/5.31.3/share/numeric/interpol.mac"
(%i4)
h:[[1,(-(13+2*b)*1)/12],[2,(-(15+2*b)*1)/7],[3,(-(17+2*b)*3)/16],[4,(-(19+2*b)*2)/9],[5,(-(21+2*b)*1)/4],[6,(-(23+2*b)*3)/11]]
(%o4)
[[1,(-2*b-13)/12],[2,(-2*b-15)/7],[3,3*(-2*b-17)/16],[4,2*(-2*b-19)/9],[5,(-2*b-21)/4],[6,3*(-2*b-23)/11]]
(%i5) ratinterpol(h,1)
(%o5)
((32*b^5+1552*b^4+30064*b^3+291128*b^2+1412514*b+2766033)*x/(32*b+16)+144)/(x^4-(2*b+53)*x^3/2+(4*b^2+128*b+1283)*x^2/4-(8*b^3+300*b^2+3974*b+19993)*x/8-(80*b^4+3440*b^3+56240*b^2+418500*b+1260063)/(16*b+8))
(%i6) ratinterpol(h,2)
(%o6) (-x^2-(2*b+11)*x/2)/(x+5)
(%i7)
h:[[1,(-1)/12],[2,(-1)/7],[3,(-3)/16],[4,(-2)/9],[5,(-1)/4],[6,(-3)/11]]
(%o7) [[1,-1/12],[2,-1/7],[3,-3/16],[4,-2/9],[5,-1/4],[6,-3/11]]
(%i8) ratinterpol(h,1)
(%o8) -x/(2*(x+5))
(%o8) "/home/user/ratinterpol-20132712.mac"
------ output with new interpol.mac ------
On 12/28/2013 03:29 AM, andre maute wrote:
> Hi list,
>
> I have the following unexpected behavior with ratinterpol.
> I would have had expected the first two ratinterpol commands
> returning something containing the unknown b.
>
> Andre
>
> ------ ratinterpol-20131227.mac ------
> display2d : false;
>
> load("interpol");
>
> h :
> [
> [1,-(2*b+13)*1/12],
> [2,-(2*b+15)*1/7],
> [3,-(2*b+17)*3/16],
> [4,-(2*b+19)*2/9],
> [5,-(2*b+21)*1/4],
> [6,-(2*b+23)*3/11]
> ];
>
> ratinterpol(h,1);
> ratinterpol(h,2);
>
> h :
> [
> [1,-1/12],
> [2,-1/7],
> [3,-3/16],
> [4,-2/9],
> [5,-1/4],
> [6,-3/11]
> ];
>
> ratinterpol(h,1);
> ------ ratinterpol-20131227.mac ------
>
> ------ output ------
> [user at localhost ~]$ maxima -b ratinterpol-20132712.mac
> Maxima 5.31.3 http://maxima.sourceforge.net
> using Lisp SBCL 1.1.8-2.fc19
> Distributed under the GNU Public License. See the file COPYING.
> Dedicated to the memory of William Schelter.
> The function bug_report() provides bug reporting information.
> STYLE-WARNING: redefining MAXIMA::$FILE_TYPE in DEFUN
> (%i1) batch("ratinterpol-20132712.mac")
>
> read and interpret file: /home/user/ratinterpol-20132712.mac
> (%i2) display2d:false
> (%o2) false
> (%i3) load("interpol")
> (%o3)
> "/home/user/opt/maxima/share/maxima/5.31.3/share/numeric/interpol.mac"
> (%i4)
> h:[[1,(-(13+2*b)*1)/12],[2,(-(15+2*b)*1)/7],[3,(-(17+2*b)*3)/16],[4,(-(19+2*b)*2)/9],[5,(-(21+2*b)*1)/4],[6,(-(23+2*b)*3)/11]]
> (%o4)
> [[1,(-2*b-13)/12],[2,(-2*b-15)/7],[3,3*(-2*b-17)/16],[4,2*(-2*b-19)/9],[5,(-2*b-21)/4],[6,3*(-2*b-23)/11]]
> (%i5) ratinterpol(h,1)
> (%o5) (1500441*x/16+144)/(x^4-55*x^3/2+1415*x^2/4-24275*x/8-579441/8)
> (%i6) ratinterpol(h,2)
> (%o6) (-x^2-13*x/2)/(x+5)
> (%i7)
> h:[[1,(-1)/12],[2,(-1)/7],[3,(-3)/16],[4,(-2)/9],[5,(-1)/4],[6,(-3)/11]]
> (%o7) [[1,-1/12],[2,-1/7],[3,-3/16],[4,-2/9],[5,-1/4],[6,-3/11]]
> (%i8) ratinterpol(h,1)
> (%o8) -x/(2*(x+5))
> (%o8) "/home/user/ratinterpol-20132712.mac"
> ------ output ------
>
> _______________________________________________
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima
>
-------------- next part --------------
/* COPYRIGHT NOTICE
Copyright (C) 2005-2012 Mario Rodriguez Riotorto
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 at
http://www.gnu.org/copyleft/gpl.html
*/
/* INTRODUCTION
This package defines some interpolation techniques.
For questions, suggestions, bugs and the like, feel free
to contact me at
mario @@@ edu DOT xunta DOT es
*/
/* Returns de input in the form of a list of pairs, ordered wrt the first */
/* coordinate. The argument must be either: */
/* a) a two column matrix, p:matrix([2,4],[5,6],[9,3]) */
/* b) a list of pairs, p: [[2,4],[5,6],[9,3]] */
/* c) a list of numbers, p: [4,6,3], in which case the abscissas will be */
/* assigned automatically to 1, 2, 3, etc. */
/* This is an uxiliary function for the 'interpol' package. */
interpol_check_input(data,funame):=
block([n,out],
if not listp(data) and not matrixp(data)
then error("Argument to '",funame,"' must be a list or matrix"),
n: length(data),
if n<2
then error("Argument to '",funame,"' has too few sample points")
elseif listp(data) and every('identity,map(lambda([x], listp(x) and length(x)=2),data))
then out: sort(data)
elseif matrixp(data) and length(data[1]) = 2
then out: sort(args(data))
elseif listp(data) and every('identity,map(lambda([x], not listp(x)),data))
then out: makelist([i,data[i]],i,1,n)
else error("Error in arguments to '",funame,"' function"),
/* controlling duplicated x's */
for i:2 thru n do
if out[i-1][1] = out[i][1]
then error("Duplicated abscissas are not allowed"),
out )$
/* Lagrangian interpolation. The argument must be either: */
/* a) a two column matrix, p:matrix([2,4],[5,6],[9,3]) */
/* b) a list of pairs, p: [[2,4],[5,6],[9,3]] */
/* c) a list of numbers, p: [4,6,3], in which case the abscissas will be */
/* assigned automatically to 1, 2, 3, etc. */
/* In cases a) and b) the pairs are ordered wrt the 1st. coordinate before any */
/* computation is made. Option: */
/* 'varname='x: the name of the independent variable */
/* Sample session: */
/* load(interpol); */
/* p:[[7,2],[8,2],[1,5],[3,2],[6,7]]$ */
/* lagrange(p); */
/* f(x):=''%; */
/* map(f,[2.3,5/7,%pi]); */
/* load(draw)$; */
/* draw2d( */
/* explicit(f(x),x,0,9), */
/* point_size = 3, */
/* points(p)) $ */
lagrange(tab,[select]) := block([n,sum:0,prod,options, defaults,ratprint:false],
tab: interpol_check_input(tab,"lagrange"),
options: ['varname],
defaults: ['x],
for i in select do(
aux: ?position(lhs(i),options),
if numberp(aux) and aux <= length(options) and aux >= 1
then defaults[aux]: rhs(i)),
if not symbolp(defaults[1])
then error("Option 'varname' is not correct"),
/* constructing the interpolating polynomial */
n: length(tab),
for i:1 thru n do(
prod: 1,
for k:1 thru n do
if k#i then prod: prod * (defaults[1]-tab[k][1]) / (tab[i][1]-tab[k][1]),
sum: sum + prod * tab[i][2] ),
sum )$
/* Characteristic function for intervals. Returns true iif */
/* z belongs to [l1, l2). This is an auxiliary function to */
/* be called from linearinterpol and cspline */
charfun2(z,l1,l2):= charfun(l1 <= z and z < l2)$
/* Linear interpolation. The argument must be either: */
/* a) a two column matrix, p:matrix([2,4],[5,6],[9,3]) */
/* b) a list of pairs, p: [[2,4],[5,6],[9,3]] */
/* c) a list of numbers, p: [4,6,3], in which case the abscissas will be */
/* assigned automatically to 1, 2, 3, etc. */
/* In cases a) and b) the pairs are ordered wrt the 1st. coordinate before any */
/* computation is made. Option: */
/* 'varname='x: the name of the independent variable */
/* Sample session: */
/* load(interpol); */
/* p:[[7,2],[8,2],[1,5],[3,2],[6,7]]$ */
/* linearinterpol(p); */
/* f(x):=''%; */
/* map(f,[2.3,5/7,%pi]); */
/* load(draw)$; */
/* draw2d( */
/* explicit(f(x),x,0,9), */
/* point_size = 3, */
/* points(p)) $ */
linearinterpol(tab,[select]) := block([n,s:0,a,b,options, defaults,ratprint:false],
tab: interpol_check_input(tab,"linearinterpol"),
options: ['varname],
defaults: ['x],
for i in select do(
aux: ?position(lhs(i),options),
if numberp(aux) and aux <= length(options) and aux >= 1
then defaults[aux]: rhs(i)),
if not symbolp(defaults[1])
then error("Option 'varname' is not correct"),
/* constructing the interpolating polynomial */
n: length(tab),
if n=2 /* case of two points */
then s: tab[2][2] + (tab[2][2]-tab[1][2]) *
(defaults[1]-tab[2][1]) /
(tab[2][1]-tab[1][1])
else for i:2 thru n do(
if i=2
then (a: 'minf, b: tab[i][1])
else if i=n
then (a: tab[i-1][1], b: 'inf)
else (a: tab[i-1][1], b: tab[i][1]),
s: s + funmake('charfun2,[defaults[1], a, b]) *
expand( tab[i][2] + (tab[i][2]-tab[i-1][2]) *
(defaults[1]-tab[i][1]) /
(tab[i][1]-tab[i-1][1]) ) ),
s )$
/* Cubic splines interpolation. The argument must be either: */
/* a) a two column matrix, p:matrix([2,4],[5,6],[9,3]) */
/* b) a list of pairs, p: [[2,4],[5,6],[9,3]] */
/* c) a list of numbers, p: [4,6,3], in which case the abscissas will be */
/* assigned automatically to 1, 2, 3, etc. */
/* In cases a) and b) the pairs are ordered wrt the 1st. coordinate before any */
/* computation is made. Options: */
/* 'd1='unknown: 1st derivative at x_1; if it is 'unknown, the second derivative */
/* at x_1 is made equal to 0 (natural cubic spline); if it is equal to a */
/* number, the second derivative is estimated based on this number */
/* 'd2='unknown: 1st derivative at x_n; if it is 'unknown, the second derivative */
/* at x_n is made equal to 0 (natural cubic spline); if it is equal to a */
/* number, the second derivative is estimated based on this number */
/* 'varname='x: the name of the independent variable */
/* Reference: this algorithm is based on 'Numerical Recipes in C', section 3.3 */
/* Sample session: */
/* load(interpol); */
/* p:[[7,2],[8,2],[1,5],[3,2],[6,7]]$ */
/* cspline(p); ==> natural cubic spline (second derivatives are zero in both extremes) */
/* f(x):=''%; */
/* map(f,[2.3,5/7,%pi]); */
/* load(draw)$; */
/* draw2d( */
/* explicit(f(x),x,0,9), */
/* point_size = 3, */
/* points(p)) $ */
/* cspline(p,d1=0,dn=0); */
/* g(x):=''%; */
/* draw2d( */
/* explicit(g(x),x,0,9), */
/* point_size = 3, */
/* points(p)) $ */
cspline(tab,[select]):= block([options, defaults, n, aux, y2, u, sig, p,
qn, un, a, b, s:0, aj, bj, cj, dj, aux,ratprint:false],
tab: interpol_check_input(tab,"cspline"),
options: ['d1, 'dn, 'varname],
defaults: ['unknown, 'unknown, 'x],
for i in select do(
aux: ?position(lhs(i),options),
if numberp(aux) and aux <= length(options) and aux >= 1
then defaults[aux]: rhs(i)),
if not numberp(defaults[1]) and defaults[1] # 'unknown
then error("Option 'd1' is not correct"),
if not numberp(defaults[2]) and defaults[2] # 'unknown
then error("Option 'dn' is not correct"),
if not symbolp(defaults[3])
then error("Option 'varname' is not correct"),
/* if tab contains only two points, linear interpolation */
n: length(tab),
if n=2 /* case of two points */
then return(ratsimp( tab[2][2] + (tab[2][2]-tab[1][2]) *
(defaults[3]-tab[2][1]) /
(tab[2][1]-tab[1][1]))),
/* constructing the interpolating polynomial */
y2: makelist(0,i,1,n),
u: makelist(0,i,1,n-1),
/* controlling the lower boundary condition */
if /*d1*/ defaults[1] = 'unknown
then (y2[1]: 0,
u[1]: 0)
else (y2[1]: -1/2,
u[1]: 3 / (tab[2][1]-tab[1][1]) *
((tab[2][2] - tab[1][2])/(tab[2][1] - tab[1][1]) - defaults[1]) ),
/* decomposition loop of the triangular algorithm */
for i:2 thru n-1 do (
sig: (tab[i][1] - tab[i-1][1]) / (tab[i+1][1] - tab[i-1][1]),
p: sig * y2[i-1] + 2,
y2[i]: (sig - 1) / p,
u[i]: (tab[i+1][2] - tab[i][2]) /(tab[i+1][1] - tab[i][1]) -
(tab[i][2] - tab[i-1][2]) /(tab[i][1] - tab[i-1][1]),
u[i]: (6 * u[i] / (tab[i+1][1] - tab[i-1][1]) - sig * u[i-1]) / p ) ,
/* controlling the upper boundary condition */
if /*dn*/ defaults[2] = 'unknown
then (qn: 0,
un: 0)
else (qn: 1/2,
un: 3 / (tab[n][1] - tab[n-1][1]) *
(defaults[2] - (tab[n][2] - tab[n-1][2]) / (tab[n][1] - tab[n-1][1]))),
y2[n]: (un - qn * u[n-1]) / (qn * y2[n-1] + 1),
/* backsubstitution loop of the tridiagonal algorithm */
for k: n-1 thru 1 step -1 do
y2[k]: y2[k] * y2[k+1] + u[k],
/* constructing the cubic splines */
for j:2 thru n do (
if j=2
then (a: 'minf, b: tab[j][1] )
else if j=n
then (a: tab[j-1][1], b: 'inf)
else (a: tab[j-1][1], b: tab[j][1]),
/* in the following sentences, defaults[3] is variable's name */
aux: (tab[j][1] - tab[j-1][1]),
aj: (tab[j][1] - defaults[3]) / aux,
bj: (defaults[3] - tab[j-1][1]) / aux,
aux: aux * aux /6,
cj: (aj^3 - aj) * aux,
dj: (bj^3 - bj) * aux,
s: s + funmake('charfun2,[defaults[3], a, b]) *
expand(aj * tab[j-1][2] + bj * tab[j][2] + cj * y2[j-1] + dj * y2[j]) ),
s )$
/* Rational interpolation, with interpolating function of the form */
/* r */
/* p x + ... + p x + p */
/* r 1 0 */
/* R(x) = ------------------------ */
/* m */
/* q x + ... + q x + q */
/* m 1 0 */
/* The 2nd. argument r is the degree of the numerator (r < sample size). The */
/* degree of the denominator is calculated as m: sample_size - r - 1. */
/* The 1st. argument must be either: */
/* a) a two column matrix, p:matrix([2,4],[5,6],[9,3]) */
/* b) a list of pairs, p: [[2,4],[5,6],[9,3]] */
/* c) a list of numbers, p: [4,6,3], in which case the abscissas will be */
/* assigned automatically to 1, 2, 3, etc. */
/* In cases a) and b) the pairs are ordered wrt the 1st. coordinate before any */
/* computation is made. Option: */
/* 'varname='x: the name of the independent variable */
/* Sample session: */
/* load(interpol)$ */
/* load(draw)$ */
/* p: [[7.2, 2.5], [8.5, 2.1], [1.6, 5.1], [3.4, 2.4], [6.7, 7.9]]$ */
/* for k:0 thru length(p)-1 do */
/* draw2d(explicit(ratinterpol(p,k),x,0,9), */
/* point_size = 3, */
/* points(p), */
/* title = concat("Degree of numerator = ",k), */
/* yrange=[0,10])$ */
ratinterpol(tab,r,[select]) :=
block([n,m,coef,unk,sol,lovtab,lov,options,defaults,ratprint:false],
tab: interpol_check_input(tab,"ratinterpol"),
options: ['varname],
defaults: ['x],
for i in select do(
aux: ?position(lhs(i),options),
if numberp(aux) and aux <= length(options) and aux >= 1
then defaults[aux]: rhs(i)),
if not symbolp(defaults[1])
then error("Option 'varname' is not correct"),
/* constructing the interpolating rational function */
n: length(tab),
if not integerp(r) or r > n-1 or r < 0
then error("Degree of numerator must be a positive integer less than sample size"),
m: n - r - 1, /* degree of denominator */
/* coef is the matrix of an homogeneous linear system */
coef: apply(matrix,
makelist(block([x,y],
[x,y]: p,
append([1],
makelist(x^k,k,1,r),
[-y],
makelist(-y*x^k,k,1,m))),
p, tab)),
unk: makelist(gensym(), k, r+m+2),
sol: map(last, linsolve(flatten(args(coef.unk)),unk)),
lovtab : listofvars(tab),
lov: listofvars(sol),
lov: listify(setdifference(setify(lov),setify(lovtab))),
sol: subst(map(lambda([z1,z2], z1=z2), lov, makelist(1, k, length(lov))), sol),
makelist(sol[k],k,1,r+1) . makelist(defaults[1]^k,k,0,r) /
makelist(sol[k],k,r+2,r+2+m) . makelist(defaults[1]^k,k,0,m) )$