axiom-math
[Top][All Lists]

Re: [Axiom-math] rational interpolation for axiom

 From: Martin Rubey Subject: Re: [Axiom-math] rational interpolation for axiom Date: Mon, 14 Jun 2004 14:03:27 +0000

root writes:
> Martin,
>
> Do you have some test cases I can add to the input subdirectory?

Not yet.

> Ideally the test cases would explain what they were trying to test
> which, so far, none of the existing src/input pamphlet files do.

Yes. I have question, too: I'd like to include the test cases in the pamphlet
file. In outher words, I'd like "document rinterp.pamphlet" to produce a third
file, namely "rinterp.input". How can I do this automatically? I want that the
documentation remains together.

> Once I've tested it I'll upload it to the world.

Please wait a little. I want to write up a third kind of interpolation that I
use frequently, namely Hermite interpolation. Furthermore, pinterp.spad is
differently structured (interfacewise) and I'm not sure which is the best. To
give you an impression, here is the second version.

Martin

\documentclass{article}
\usepackage{axiom,amsmath,url}
\begin{document}
\author{Martin Rubey}
\maketitle
\begin{abstract}
Rational Interpolation
\end{abstract}
\section{Introduction}
This file contains a crude na\"ive implementation of rational interpolation,
where the coefficients of the rational function are in any given field.

\section{Questions and Outlook}
\begin{itemize}
\item Maybe this file should be joined with pinterp.spad, where polynomial
Lagrange interpolation is implemented. This version parallels the structure of
\item There are probably better ways to implement rational interpolation. Maybe
\url{http://www.cs.ucsb.edu/~omer/personal/abstracts/rational.html} contains
something useful, but I don't know.
\end{itemize}

<<*>>=
<<RINTERP Exports>>
<<RINTERP Implementation>>
@

\section{RationalInterpolation}

)abbrev package RINTERPA RationalInterpolationAlgorithms
++ Description:
++ This package exports rational interpolation algorithms
RationalInterpolationAlgorithms(F, P): Cat == Body   where
F: Field
P: UnivariatePolynomialCategory(F)
@

<<RINTERP Exports>>=
Cat == with
RationalInterpolation: (List F, List F, NonNegativeInteger,
NonNegativeInteger) -> Fraction P
@

The implementation sets up a system of linear equations and solves it.
<<RINTERP Implementation>>=
RationalInterpolation(xlist, ylist, m, k) ==
@

First we check whether we have the right number of points and values. Clearly
the number of points and the number of values must be identical. Note that we
want to determine the numerator and denominator polynomials only up to a
factor. Thus, we want to determine $m+k+1$ coefficients, where $m$ is the degree
of the polynomial in the numerator and $k$ is the degree of the polynomial in
the denominator.

In fact, we could also leave -- for example -- $k$ unspecified and determine it
as $k=[[#xlist]]-m-1$: I don't know whether this would be better.
<<RINTERP Implementation>>=
#xlist ^= #ylist =>
error "Different number of points and values."
#xlist ^= m+k+1 =>
error "wrong number of points"
@

The next step is to set up the matrix. Suppose that our numerator polynomial is
$p(x)=a_0+a_1x+\dots+a_mx^m$ and that our denominator polynomial is
$q(x)=b_0+b_1x+\dots+b_mx^m$. Then we have the following equations, writing $n$
for $m+k+1$:

\begin{align*}
p(x_1)-y_1q(x_1)&=a_0+a_1x_1+\dots +a_mx_1^m-y_1(b_0+b_1x_1+\dots
+b_kx_1^k)=0\\
p(x_2)-y_2q(x_2)&=a_0+a_1x_2+\dots +a_mx_2^m-y_2(b_0+b_1x_2+\dots
+b_kx_2^k)=0\\
&\;\;\vdots\\
p(x_n)-y_nq(x_n)&=a_0+a_1x_n+\dots +a_mx_n^m-y_n(b_0+b_1x_n+\dots +b_kx_n^k)=0
\end{align*}

This can be written as
$$\begin{pmatrix} 1&x_1&\dots&x_1^m&-y_1&-y_1x_1&\dots&-y_1x_1^k\\ 1&x_2&\dots&x_2^m&-y_2&-y_2x_2&\dots&-y_2x_2^k\\ \vdots\\ 1&x_n&\dots&x_n^m&-y_n&-y_nx_n&\dots&-y_nx_2^k \end{pmatrix} \begin{pmatrix} a_0\\a_1\\\vdots\\a_m\\b_0\\b_1\\\vdots\\b_k \end{pmatrix}=\mathbf 0$$
We generate this matrix columnwise:
<<RINTERP Implementation>>=
tempvec: List F := [1 for i in 1..(m+k+1)]

collist: List List F := cons(tempvec,
[(tempvec := [tempvec.i * xlist.i _
for i in 1..(m+k+1)]) _
for j in 1..max(m,k)])

collist := append([collist.j for j in 1..(m+1)], _
[[- collist.j.i * ylist.i for i in 1..(m+k+1)] _
for j in 1..(k+1)])
@
Now we can solve the system:
<<RINTERP Implementation>>=
res: List Vector F := nullSpace((transpose matrix collist) _
::Matrix F)
@

Note that it may happen that the system has several solutions. In this case,
some of the data points may not be interpolated correctly. However, the
solution is often still useful, thus we do not signal an error.

<<RINTERP Implementation>>=
if #res~=1 then output("Warning: unattainable points!" _
::OutputForm)$OutputPackage @ In this situation, all the solutions will be equivalent, thus we can always simply take the first one: <<RINTERP Implementation>>= reslist: List List P := _ [[monomial((res.1).(i+1),i) for i in 0..m], _ [monomial((res.1).(i+m+2),i) for i in 0..k]] @ Finally, we generate the rational function: <<RINTERP Implementation>>= reduce((_+),reslist.1)/reduce((_+),reslist.2) <<*>>= )abbrev package RINTERP RationalInterpolation ++ Description: ++ This package exports interpolation algorithms RationalInterpolation(xx, F): Cat == Body where xx: Symbol F: Field UP ==> UnivariatePolynomial SUP ==> SparseUnivariatePolynomial Cat == with interpolate: (Fraction UP(xx,F), List F, List F, _ NonNegativeInteger, NonNegativeInteger) _ -> Fraction UP(xx,F) interpolate: (List F, List F, NonNegativeInteger, NonNegativeInteger) _ -> Fraction SUP F Body == add RIA ==> RationalInterpolationAlgorithms interpolate(qx, lx, ly, m, k) == px := RationalInterpolation(lx, ly, m, k)$RIA(F, UP(xx, F))
elt(px, qx)

interpolate(lx, ly, m, k) ==
RationalInterpolation(lx, ly, m, k)\$RIA(F, SUP F)
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}