axiom-developer
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

## [Axiom-developer] [RationalInterpolation]

 From: kratt6 Subject: [Axiom-developer] [RationalInterpolation] Date: Tue, 22 Mar 2005 04:44:32 -0600

Changes
http://page.axiom-developer.org/zope/mathaction/RationalInterpolation/diff
--

??changed:
-The package below implements rational interpolation.
-
-It requires the following previously compiled package:
Introduction

This file contains an implementation of rational interpolation, where the data
points are element of any integral domain.

Questions and Outlook

- Maybe this file should be joined with pinterp.spad, where polynomial
Lagrange interpolation is implemented. This version parallels the structure
of pinterp.spad closely. This also answers comments and questions from
wyscc. He remarked

- Abbreviations for a constructor should be limited to 7 letters (not 8).
The system occasionally adds the 8th character to a package for internal
use.

- Function names begin with a lower case, so RationalInterpolation should
be rationalInterpolation, or better, rationalInterpolate.

- Regarding the types I used for the values, wyscc remarked

- If we are doing a rational interpolation, presumably the values are
rational, so it does not make sense to require the $y$-coordinates of
inputs be integral. On the other hand, as in the above example, if one
uses 'FRAC INT', problems can arise when this package is combined with
other packages that constructs the quotient field of the parameter domain
'F' because Axiom does not like constructing 'FRAC FRAC INT'.

Note however, that the package would rather construct the type 'FRAC SUP
FRAC INT', so this problem should not occur. Moreover, there are situations
- for example in the package [Guess], where we want to interpolate values
from an IntegralDomain. Of course we could first convert them to the
quotient field, however, the current approach seems more natural to me.

- Finally, wyscc asked:
If <code>p(xx) = interpolate(lx, ly, m, k)</code>, what is the purpose of
<code>elt(px, qx) = p(qx)</code>, the composition of <code>p(xx)</code> and
<code>qx</code>, especially when <code>qx</code> is from <code>FRAC UP(xx,
F)</code> instead of from just <code>F</code>? and why is this function
(the composition) also called <code>interpolate</code>?

I do not really know - apart from a very superficial level: Clearly, the
second function was intended to let the user easily plug in values into the
interpolated function. I don't find this sensible and I would be happy to
change it. Indeed, this would also get rid of the first parameter to
'RINTERP', which is quite a nuisance.

I think we should agree on a general interface for interpolation
algorithms, and mark 'PINTERP' as obsolete. By the way, it seems that
'RINTERP' is faster, too.

- There are probably better ways to implement rational interpolation. Maybe
http://www.cs.ucsb.edu/~omer/personal/abstracts/rational.html
contains something useful, but I don't know.

- For those who speak german,
http://www.num.math.uni-goettingen.de/schaback/teaching/numath.ps
contains quite a bit of information.

- This implementation of rational interpolation neither takes care of
unattainable points, nor does it check whether the values of the
$x$-coordinates are all distinct.

- Comments welcome!

??changed:
-)lib RINTERPA
)abbrev package RINTERPA RationalInterpolationAlgorithms
++ Description:
++ This package exports rational interpolation algorithms
RationalInterpolationAlgorithms(F, P): Cat == Body   where
F: IntegralDomain
P: UnivariatePolynomialCategory(F)
Cat == with
RationalInterpolation: (List F, List F, NonNegativeInteger,
NonNegativeInteger)
-> Fraction P
+++ We assume that the elements of the first list are all distinct.
+++ If they are not, division by zero might occur.

Body == add
RationalInterpolation(xlist, ylist, m, k) ==
#xlist ^= #ylist =>
error "Different number of points and values."
#xlist ^= m+k+1 =>
error "wrong number of points"
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)])
resspace: List Vector F := nullSpace((transpose matrix collist) _
::Matrix F)
reslist: List List P := _
[[monomial((resspace.1).(i+1), i) for i in 0..m], _
[monomial((resspace.1).(i+m+2), i) for i in 0..k]]

reduce((_+), reslist.1)/reduce((_+), reslist.2)

??changed:
-Next PolynomialCommonDenominator
-
-Example (added by wyscc):

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.

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{eqnarray*}
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{eqnarray*}

This can be written as
\begin{equation*}
\begin{bmatrix}
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{bmatrix}
\begin{bmatrix}
a_0\\a_1\\\vdots\\a_m\\b_0\\b_1\\\vdots\\b_k
\end{bmatrix}=\mathbf 0
\end{equation*}

We generate this matrix columnwise, then we can solve the system using
'nullSpace'.

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.

Since all the solutions of 'nullSpace' will be equivalent, we can always
simply take the first one. Finally, we return the rational function.

Examples

To conclude we present some examples. To begin with, the following
interpolation
illustrates the concept of unattainable points:

??changed:
-f(x)== (x^3+5*x-3)/(x^2-3)
-f(x)
-xlist:List FRAC INT :=[1/2, 4, 1/6, 8, 1/10, 12]
-ylist :=[f(x) for x in xlist]
interpolate([q,q^2,q^3],[0,x^1,x^2],0,2)$RINTERP(qn, FRAC POLY INT) ++added: \begin{axiom} f(x) == (x^3+5*x-3)/(x^2-3) xlist := [1/2, 4, 1/6, 8, 1/10, 12] ylist := [f(x) for x in xlist] interpolate(xlist, ylist, 3, 2)$RINTERP('x, FRAC INT)
interpolate(1/6::FRAC UP(x,FRAC INT), xlist, ylist, 3, 2)$RINTERP('x,FRAC INT) \end{axiom} ??changed: -dom:=DMP([z],INT); -g:FRAC dom -> FRAC dom dom := DMP([z],INT); g: FRAC dom -> FRAC dom; ??changed: -xxlist:List FRAC dom:=[1/(2*z), 4*z, 1/(6*z), 8*z, 1/(10*z), 12*z] -yylist:=[g(x) for x in xxlist] ---RationalInterpolation(xxlist, yylist, 3::NNI, 2::NNI)$RINTERPA(FRAC dom,
UP(x, FRAC dom))
-\end{axiom}
-
-\begin{axiom}
-interpolate(xlist, ylist, 3, 2)$RINTERP('x, FRAC INT) -interpolate(1/6::FRAC UP(x,FRAC INT), xlist, ylist, 3,2)$RINTERP('x,FRAC INT)
xxlist: List FRAC dom := [1/(2*z), 4*z, 1/(6*z), 8*z, 1/(10*z), 12*z]
yylist := [g(x) for x in xxlist]

--removed:
-
-Question: If <code>p(xx) = interpolate(lx, ly, m, k)</code>, what is the
purpose of
-<code>elt(px, qx) = p(qx)</code>, the composition of <code>p(xx)</code> and
<code>qx</code>, especially when <code>qx</code> is from <code>FRAC UP(xx,
F)</code> instead of from just <code>F</code>? and why is this function (the
composition) also called <code>interpolate</code>?
-

--
forwarded from http://page.axiom-developer.org/zope/mathaction/address@hidden



reply via email to

 [Prev in Thread] Current Thread [Next in Thread]