Re: [Axiom-mail] Another beginner's question: Lagrange interpolation
Dear Alasdair, Dear Vadim!
Vadim: sorry I haven't sent you the paper by Gemignani yet. It turned out that
I do not have a pdf. If you are still interested, please send me a postal
adress and I send you a collection of papers today. See below on a tiny
demonstration what FFFG can do.
Alasdair McAndrew [EMAIL PROTECTED] writes:
Just to explore Axiom a little further, I was trying to do a Lagrange
interpolation on two lists:
LagrangeInterpolation([1,2,3],[4,5,-6])
This doesn't work - I get an error message which I (as a beginner) don't
understand. So then:
1) How do I do this?
see Bill's mail or below.
2) Without asking this mail group, where can I find this information in a
manner which I can understand?
If you do not want to ask on axiom-mail, you should ask on axiom-math. :-)
Seriously, I think the best and most rapid source of information is
axiom-math or axiom-mail for how do I do interpolation, symmetric functions,
limits, plot graphs, etc.
axiom-developer for why doesn't axiom compile on my machine and I have a
week spare time, and I'm sooo bored, what shall I do.
I'd advise you to spend roughly fifteen minutes
* browsing HyperDoc,
* using )what things interpol
* searching the pdf version of the book
* searching the mail archives and the wiki
If you are not successful, send mail.
-- a more scientific answer -
interpolation:
Lagrange interpolation as implemented in Axiom is fine for toy demonstrations.
But if you are interested in interpolation in general, I suggest that you use
my implementation of Bernhard Beckermann and George Labahn's FFFG (Fast
Fraction Free Gaussian) elimination algorithm. It is faster in exact domains
and much more general. The same file also features (in a separate package)
polynomial interpolation using Newton, albeit in an extremely specialised
setting:
++ \spad{newton}(l) returns the interpolating polynomial for the values
++ l, where the x-coordinates are assumed to be [1,2,3,...,n] and the
++ coefficients of the interpolating polynomial are known to be in the
++ domain F. I.e., it is a very streamlined version for a special case of
++ interpolation.
You will have to install my package manually though, or use wh-sandbox. Here
are the docstrings:
interpolate: (List Fraction D, List Fraction D, NonNegativeInteger)
- Fraction SUP D
++ \spad{interpolate(xlist, ylist, deg} returns the rational function with
++ numerator degree \spad{deg} that interpolates the given points using
++ fraction free arithmetic.
generalInterpolation: (List D, CoeffAction,
Vector V, List NonNegativeInteger) - Matrix SUP D
++ \spad{generalInterpolation(C, CA, f, eta)} performs Hermite-Pade
++ approximation using the given action CA of polynomials on the elements
++ of f. The result is guaranteed to be correct up to order
++ |eta|-1. Given that eta is a normal point, the degrees on the
++ diagonal are given by eta. The degrees of column i are in this case
++ eta + e.i - [1,1,...,1], where the degree of zero is -1.
++
++ The first argument C is the list of coefficients c_{k,k} in the
++ expansion x^k z g(x) = sum_{i=0}^k c_{k,i} x^i g(x).
++
++ The second argument, CA(k, l, f), should return the coefficient of x^k
++ in z^l f(x).
So, polynomial interpolation is interpolate([1,2,3],[4,5,-6], 3), Hermite-Pade
is, for example
(1) - f :=reduce(+, [monomial(random 40, i)$SUP INT for i in 0..5])
5 4 3 2
(1) 12? + 25? + 29? + 13? + 20? + 1
Type: SparseUnivariatePolynomial Integer
(2) - F == FFFG(INT, SUP INT)
Type: Void
(3) - eta := [2,1,1]
(3) [2,1,1]
Type: List PositiveInteger
(4) - n := reduce(+, eta)
(4) 4
Type: PositiveInteger
(5) - M := generalInterpolation(DiffC(n)$F, DiffAction$F, [f, D f, D D f],
eta)$F
(5)
+ 2 +
|- 1045836? + 939076? + 1012588 3237618? + 3755082 16407212? + 9862196 |
| |
|133656 - 1045836? + 2994301116948 |
| |
+ - 141758 - 374757- 1045836? - 1238506+
Type: Matrix SparseUnivariatePolynomial Integer
(6) - [f, D f, D D f] * M
(6)
76 54
[- 12550032? - 14876988? + 5298712? + 46971396? ,
65
RE: [Axiom-mail] Another beginner's question: Lagrange interpolation
On April 13, 2007 10:07 AM Martin Rubey wrote: Alasdair McAndrew writes: I wish I could get more of a handle on Axiom more quickly, but its behaviour still seems tantalizingly just out of reach. I felt like that for about the first year or so ... :-) would be interesting if you could be more precise here... I find that the interpreter is quite clever about finding appropriate types. But of course, often it needs help. And certainly, a (conservative) redesign of the hierarchy would be in order, but this also depends on the future of Aldor.. ... However, for me at least, the problem of getting the types/domains and so on correct for the use of some things (like LagrangeInterpolation) seems unnecessarily difficult. LagrangeInterpolation is really very poorly done, especially with respect to the usage of domains. I am not entirely sure I agree with Martin that this code is very poorly done. Perhaps that is true, from the point of view of the Axiom interpreter because it requires that the user provide some additional information - the type for the underlying field. But given that, the interpreter is pretty smart about deciding on the types. It is actually much easier than my previous reply might have suggested: (1) - )expose PolynomialInterpolationAlgorithms PolynomialInterpolationAlgorithms is now explicitly exposed in frame initial (2) - LagrangeInterpolation([1,2,3],[4,5,-6])@SUP(FRAC INT) 2 (3) - 6? + 19? - 9 Type: SparseUnivariatePolynomial Fraction Integer The critical part the interpreter was missing was just the choice of field. But without a good knowledge of the Axiom type system this solution would probably be rather hard for a beginner to anticipate. ... And of course I want to have a go at some programming... I've just written (the first version of) a z-transforms package for Maxima, and it seems to me that such a thing should be relatively straightforward in Axiom, using its rules mechanism. Then there'd be a method for solving difference equations! I would like to give you an unsolicited caution regarding the use of pattern matching rules in Axiom. Although the application of rules work fairly well within the Axiom 'Expression' domain really this sort of symbolic manipulation of expressions is rather foreign and unnatural in Axiom. People who program in Python have invented a name of the style of programming that is natural in that language, i.e. pythonic (sic). In Axiom we can identify an analogous style that we might call Axiomic or maybe Axiomatic. :-) The point being that most things in Axiom are done in an object- oriented, algebraic manner -- and not by the simple manipulation of expressions. As a long time user of Maple, Maxima and other symbol-oriented systems, I think the hardest part about learning to use Axiom for me was to understand this change in paradigm. Most beginning users of Axiom think of types as just getting in the way, but in the end that is really where the strength of Axiom actually lies. One thing that the object-orientation buys you (and maybe costs you) is the importance of understanding the existing Axiom library (which currently contains nearly 1,300 different abstract mathematical and data types. To get the greatest advantage and do the least amount of programming it is essential that you think about how what you intend to write fits into the existing mathematical hierarchy. This means more time spent up front analysing and generalizing your algorithm but ultimately less time actually programming. In your case perhaps you need to think about how the z-transform relates to other integral transforms that already exist in Axiom. For example: (5) - )wh package transform Packages --- Packages with names matching patterns: transform INVLAPLA InverseLaplaceTransform LAPLACE LaplaceTransform (5) - )sh LAPLACE LaplaceTransform(R: Join(EuclideanDomain,OrderedSet,CharacteristicZero, RetractableTo Integer,LinearlyExplicitRingOver Integer), F: Join(TranscendentalFunctionCategory,PrimitiveFunctionCategory, AlgebraicallyClosedFunctionSpace R)) is a package constructor Abbreviation for LaplaceTransform is LAPLACE This constructor is exposed in this frame. Issue )edit c:/Program Files/axiom/mnt/windows/../../src/algebra/LAPLACE.spad to see algebra source code for LAPLACE --- Operations laplace : (F,Symbol,Symbol) - F (5) - laplace(x+1,x,y) y + 1 (5) - 2 y Type: Expression Integer -- Of course understanding what goes into the definition of the LaplaceTransform package: EuclideanDomain OrderedSet CharacteristicZero, RetractableTo Integer LinearlyExplicitRingOver Integer TranscendentalFunctionCategory PrimitiveFunctionCategory AlgebraicallyClosedFunctionSpace
Re: [Axiom-mail] Another beginner's question: Lagrange interpolation
Dear Bill, Alasdair, Bill Page [EMAIL PROTECTED] writes: I would like to give you an unsolicited caution regarding the use of pattern matching rules in Axiom. Although the application of rules work fairly well within the Axiom 'Expression' domain really this sort of symbolic manipulation of expressions is rather foreign and unnatural in Axiom. Although pattern matching and rule-based algorithms *may* not be the best way to program for axiom (and I'm not even sure about that), I believe meanwhile that writing some programs in the scripting language will give you a very fast start into axiom. Martin ___ Axiom-mail mailing list [EMAIL PROTECTED] http://lists.nongnu.org/mailman/listinfo/axiom-mail
RE: [Axiom-mail] Another beginner's question: Lagrange interpolation
Alasdair, Thanks for the question! On Sent: April 12, 2007 9:45 PM Alasdair McAndrew wrote: Just to explore Axiom a little further, I was trying to do a Lagrange interpolation on two lists: LagrangeInterpolation([1,2,3],[4,5,-6]) This doesn't work - I get an error message which I (as a beginner) don't understand. This is what I get: (1) - LagrangeInterpolation([1,2,3],[4,5,-6]) There are no exposed library operations named LagrangeInterpolation but there is one unexposed operation with that name. Use HyperDoc Browse or issue )display op LagrangeInterpolation to learn more about the available operation. Cannot find a definition or applicable library operation named LagrangeInterpolation with argument type(s) List PositiveInteger List Integer Perhaps you should use @ to indicate the required return type, or $ to specify which version of the function you need. So then: 1) How do I do this? I recommend first of all using the suggested command: (1) - )display op LagrangeInterpolation There is one unexposed function called LagrangeInterpolation : [1] (List D3,List D3) - D1 from PolynomialInterpolationAlgorithms( D3,D1) if D3 has FIELD and D1 has UPOLYC D3 2) Without asking this mail group, where can I find this information in a manner which I can understand? Hmmm what should I assume about what you can understand? That is a difficult problem and the answer will change over time. With the above information you should be able to see that at least that: - the LagrangeInterpolation function does exist in Axiom but it is not exposed - that it expects two lists as arguments - and these lists must consist of elements of a Field - and that it will return a polynomial over this Field - and finally that what you passed as arguments to this function did not meet these requirements. The first thing I would try is to look in the Axiom book by using the search function in your PDF reader. No luck. But at least you should be able to read there about what the phrase not exposed means. Essentially it means that you must call the function explicitly from the package where it is defined. The second thing you might try is to look this function up in HyperDoc using Browse, enter the name LagrangeInterpolation and click Operations. Unfortunately this only tells you what you know already plus the fact that it is not documented yet... You should also ask Axiom about 'PolynomialInterpolationAlgorithms' (1) - )show PolynomialInterpolationAlgorithms PolynomialInterpolationAlgorithms(F: Field,P: UnivariatePolynomialCategory F) is a package constructor Abbreviation for PolynomialInterpolationAlgorithms is PINTERPA This constructor is not exposed in this frame. Issue )edit c:/Program Files/axiom/mnt/windows/../../src/algebra/PINTERPA.spad to see algebra source code for PINTERPA --- Operations LagrangeInterpolation : (List F,List F) - P So now take another look at the last point and what you actually passed to the function: List PositiveInteger List Integer So Axiom thought the first list [1,2,3] was a list of PositiveInteger and the second list [4,5,-6] was a list of Integer. Lets ask Axiom if these are Fields: (1) - PositiveInteger has Field (1) false Type: Boolean (2) - Integer has Field (2) false Type: Boolean but at least (3) - Fraction Integer has Field (3) true Type: Boolean Oh! So how can I tell Axiom that these numbers are actually from some Field? Try this: (4) - [1,2,3]::List Fraction Integer (4) [1,2,3] Type: List Fraction Integer Ok, now lets try LagrangeInterpolation again: (9) - A:=[1,2,3]::List Fraction Integer (9) [1,2,3] Type: List Fraction Integer (10) - B:=[4,5,-6]::List Fraction Integer (10) [4,5,- 6] Type: List Fraction Integer But how does one call the specific package where this function is located? The command ')show PolynomialInterpolationAlgorithms' gave use this information: PolynomialInterpolationAlgorithms(F: Field,P: UnivariatePolynomialCategory F) So now in my example above the Field is 'Fraction Integer' but it may not be entirely obvious to you what UnivariatePolynomialCategory F means. To make a long story short: Axiom groups it's mathematical domains (types) into groups called Category. This is almost (but not quite) like a category in mathematics. Anyway let me suggest one such domain that will work for us: The univariate polynomials over the variable x with coefficients in the same field:
