diff options
| author | dos-reis <gdr@axiomatics.org> | 2007-08-14 05:14:52 +0000 |
|---|---|---|
| committer | dos-reis <gdr@axiomatics.org> | 2007-08-14 05:14:52 +0000 |
| commit | ab8cc85adde879fb963c94d15675783f2cf4b183 (patch) | |
| tree | c202482327f474583b750b2c45dedfc4e4312b1d /src/algebra/rinterp.spad.pamphlet | |
| download | open-axiom-ab8cc85adde879fb963c94d15675783f2cf4b183.tar.gz | |
Initial population.
Diffstat (limited to 'src/algebra/rinterp.spad.pamphlet')
| -rw-r--r-- | src/algebra/rinterp.spad.pamphlet | 150 |
1 files changed, 150 insertions, 0 deletions
diff --git a/src/algebra/rinterp.spad.pamphlet b/src/algebra/rinterp.spad.pamphlet new file mode 100644 index 00000000..41fa6c59 --- /dev/null +++ b/src/algebra/rinterp.spad.pamphlet @@ -0,0 +1,150 @@ +\documentclass{article} +\usepackage{axiom} +\begin{document} +\title{\$SPAD/src/algebra rinterp.spad} +\author{Martin Rubey} +\maketitle +\begin{abstract} +Rational Interpolation +\end{abstract} +\eject +\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. I have a second version that parallels + the structure of pinterp.spad closely. +\item 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. +\item Comments welcome! +\end{itemize} + +\section{RationalInterpolation} + +<<RINTERP Header>>= +)abbrev package RINTERP RationalInterpolation +++ Description: +++ This package exports rational interpolation algorithms +RationalInterpolation(xx,F): Exports == Implementation where + xx: Symbol + F: Field +@ + +<<RINTERP Exports>>= + Exports == with + interpolate: (List F, List F, NonNegativeInteger, + NonNegativeInteger) -> Fraction Polynomial F +@ + +The implementation sets up a system of linear equations and solves it. +<<RINTERP Implementation>>= + Implementation == add + interpolate(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$: +\noindent +$$ +\begin{array}{rl} + 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{array} +$$ +This can be written as +$$ +\left[ +\begin{array}{cccccccc} +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{array} +\right] +\left[ +\begin{array}{c} +a_0\\a_1\\\vdots\\a_m\\b_0\\b_1\\\vdots\\b_k +\end{array} +\right] +=\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 Polynomial F := _ + [[(res.1).(i+1)*(xx::Polynomial F)**i for i in 0..m], _ + [(res.1).(i+m+2)*(xx::Polynomial F)**i for i in 0..k]] +@ +Finally, we generate the rational function: +<<RINTERP Implementation>>= + reduce((_+),reslist.1)/reduce((_+),reslist.2) +@ +\section{Rational Interpolation Code} +<<package RINTERP RationalInterpolation>>= +<<RINTERP Header>> +<<RINTERP Exports>> +<<RINTERP Implementation>> +@ +<<*>>= +<<RINTERP Header>> +<<RINTERP Exports>> +<<RINTERP Implementation>> +@ +\eject +\begin{thebibliography}{99} +\bibitem{1} nothing +\end{thebibliography} +\end{document} |