Dear Barton,
I have been very glad to hear that you are updating your so interesting
and useful specfun package for special functions including all classical
systems of orthogonal polymomials. Please, accept my sincere
congratulations for this effort, particularly for your effort to include
running error analysis in the numerical values of the functions. This
will be extremely helpful especially inside a computer algebra system
as is the case with Maxima.
With this message I will try to bring to your attention my views with
respect to the possibilities of new features in this already excellent
package of yours in reply to your related request. Please, accept my
sincere apologies for the delay. (My duties do not permit me to be
more active with respect to Maxima.) Please, be also so kind to
forgive me for any incorrect, inappropriate or just irrelevant
comments below. Many thanks in advance!
> Barton Willis wrote:
> I'm seeking input on two problems and taking suggestions for new
> features.
> If you have suggestions for the new version, let me know.
Naturally, I cannot help with respect to the two problems you
mention although other colleagues have been already interested
and helped (at least for the first of them). My suggestions for the
further improvement of your already excellent package concern:
1. At first, I would like (very much indeed) that the ratsimp command
could be used in most functional evaluations. Please, compare
legendre_p(10,x);
legendre_p(10,x), ratsimp;
legendre_q(10,x);
legendre_q(10,x), ratsimp;
or
chebyshev_t(10,x);
chebyshev_t(10,x), ratsimp;
To be sincere, so far I am introducing the ratsimp command myself
in the Maxima computations, but this should be made automatically.
(This is my opinion and I may be wrong in it. I am still a novice in
using Maxima!)
2. Another point which I feel is not included now is the derivation of
series of functions based on orthogonal poynomials, e.g. Chebyshev
or Legendre series. Maple has the chebyshev command. I do not
know what is present in Maxima now. (I have to study!) In any case,
I suggest this addition to the package for the classical series
(approximate series too) especially for the main types of Chebyshev
series (if not already present, I doubt.) I mean both exact series and
series based on numerical integration.
3. Another (obvious) possibility is the addition of several more special
functions. I do not have a concrete view about which special functions
would be best for Maxima. Perhaps the exponential integral, the
sine and cosine integrals and so many additional functions. (The good
thing would be also that these functions can appear in the integration
algorithms as is already the case with the error function erf, but this
is
a difficult task.)
4. As far as I have observed, beyond the classical orthogonal
polynomials
there is also present the Legendre function of the second kind Q[n](x)
with x in (-1, 1). I do not know if it would be interesting to include
the
whole set of functions of the second kind for the classical systems of
orthogonal polynomials already present in your specfun package both
for x in [a, b] (the orthogonality interval itself, as is already the
case with
the Chebyhev and the Legendre polynomials) and outside it in the=20
complex plane z = x+iy too. This is a natural extension of your
package.
A function of the second kind may be defined as
integrate(orthogonal_polynomial(t)/(t-z), t, a, b)
with [a, b] the orthogonality interval. When z is in (a, b), we have a
Cauchy-type principal value integral. (This is already essentially the
case with Q[n](x).)
5. Perhaps, it could also be possible for you to include the algorithm
for
the construction of the system of orthogonal polynomials corresponding=20
to an arbitrary positive weight function w(x) on an interval [a, b]. (In
fact,
there are so many possibilities for extensions to your specfun package!)
6. In the past, before thirty to twenty years I used the classical
Gaussian
integration rules corresponding to all systems of orthogonal polynomials
in your package by using Fortran. Therefore, I would be somewhat
interested to see this possibility in Maxima too (if not already
present,
again I am not sure). Having your specfun pakage available, this is now
a rather easy task. For example, for the classical Gauss-Legendre rule
a first (and so incomplete!) attempt might be the following one:
PRELIMINARY COMMANDS
(C1) display2d : false$
(C2) load(specfun)$
THE AUXILIARY FUNCTION FOR THE NODES AND WEIGHTS
(GAUSS-LEGENDRE CASE ONLY)
(C3) gauss_legendre_nodes_weights(n) := block([x, p, dp, r],
p : ratsimp(legendre_p(n,x)),
dp : ratsimp(legendre_p(n+1,x)),
r : bfloat(realroots(p, 10^(-16))),
nodes : makelist(bfloat(rhs(r[k])), k,1,n),
weights :
makelist(2*(1-nodes[k]^2)/((n+1)*subst(x=nodes[k],dp))^2, k,1,n),
"The nodes and weights of the Gauss-Legendre rule found"
)$
EXAMPLE FOR THE NODES AND WEIGHTS (with n = 10)
(C4) gauss_legendre_nodes_weights(10);
(D4) "The nodes and weights of the Gauss-Legendre rule found"
(C5) nodes;
(D5) [-9.739065285171717B-1,-8.650633666889845B-1,-6.794095682990244B-1,
-4.333953941292472B-1,-1.488743389816312B-1,1.488743389816312B-1,
4.333953941292472B-1,6.794095682990244B-1,8.650633666889845B-1,
9.739065285171717B-1]
(C6) weights;
(D6) [6.667134430869263B-2,1.494513491505657B-1,2.190863625159814B-1,
2.69266719309996B-1,2.955242247147529B-1,2.955242247147529B-1,
2.69266719309996B-1,2.190863625159814B-1,1.494513491505657B-1,
6.667134430869263B-2]
TESTS (for m = 20 > 2n-1 = 19 the quadrature rule is not exact any
more!)
(C7) tests : makelist(["m =",m, integrate(x^m, x,-1,1) -
sum(weights[k]*nodes[k]^m, k,1,10)], m,1,20);
(D7) [["m =",1,5.204170427930421B-18], ["m
=",2,1.451616604697392B-14],
["m =",3,1.734723475976807B-18], ["m
=",4,8.916478666520788B-15],
["m =",5,8.673617379884035B-19], ["m
=",6,4.961309141293668B-15],
["m =",7,8.673617379884035B-19], ["m
=",8,2.144118216307334B-15],
["m =",9,8.673617379884035B-19], ["m
=",10,1.318389841742373B-16],
["m =",11,-1.734723475976807B-18], ["m
=",12,-1.283695372222837B-15],
["m =",13,0.0B0],["m =",14,-2.275957200481571B-15],
["m =",15,-8.673617379884035B-19], ["m
=",16,-2.945560462208618B-15],
["m =",17,0.0B0],["m =",18,-3.382710778154774B-15], ["m
=",19,0.0B0],
["m =",20,2.925590327091873B-6]]
THE MAIN FUNCTION
(C8) gauss_legendre(function, n) :=
sum(weights[k]*subst(x=nodes[k], function), k,1,n)$
EXAMPLES OF THE COSTRUCTED RULE
(C9) bfloat(integrate(1/(1+x^2), x,-1,1)) =
gauss_legendre(1/(1+x^2),10);
(D9) 1.570796326794897B0 = 1.570796270223256B0
(C10) bfloat(integrate(1/(10+x^2), x,-1,1)) =
gauss_legendre(1/(10+x^2),10);
(D10) 1.937068164680779B-1 = 1.937068164680757B-1
(C11) bfloat(integrate(exp(x), x,-1,1)) = gauss_legendre(exp(x),10);
(D11) 2.350402387287603B0 = 2.350402387287572B0
I believe that in this preliminary version the results are more or less
satisfactory. In case you may be interested in Gauss (and Lobatto)
rules for your specfun package (all systems of orthogonal polynomials
in it), I would be very glad to prepare and send you the analogous
Maxima code for the additional rules too so that you can proceed. In
case you have some serious reason not to be interested, possibly I
might consider the possibility of preparing a related separate package
myself (if I can find the related time some time in the future,
sincerely
I doubt) although my experience with Maxima is presently rather
insufficient.
(Naturally, by no means can I exclude that Gaussian rules are already
present in some package of Maxima!)
> I think orthopoly is better name than is specfun; would it be okay
> to change it?
I think this change is not very appropriate at first since it is a
change
in a package already well known and, moreover, since your specfun
package already includes functions which are not orthogonal polynomials.
In Maple there is (or, at least, there was) an orthopoly package, but
this package concerned only the classical orthogonal polynomials
and not additional functions (such as spherical Hankel and Bessel
functions). I could possibly agree with the functions of the second
kind I mentioned previously, but not for additional functions. There-
fore, I feel the name can remain the same.
Of course, if in the future you wish to devote the package to orthogonal
polynomials, perhaps it could be recommended to change its name
to orthopoly, reserving the name specfun to a package devoted to
non-polynomial special functions. Naturally, as Ray already mentioned,
the choice is yours!
Finally, one quite minor point: The phrase in your specfun package
that the Jacobi polynomials "When a <= -1 or b <= -1, the
polynomials
are defined, but they don't form an orthogonal set." is clearly CORRECT
with respect to the classical meaning of integrals, but if one proceeds
to accept the finite-part sense of an integral, the Jacobi polynomials
still form an orthogonal set (for any values of a and b) and,
consequently,
the related Gauss-Jacobi quadrature rule can also be used (even for
a <= -1 and/or b <= -1).
In fact, Maxima accepts Cauchy-type principal value integrals. Example:
(C12) integrate(1/x, x,-1,1);
Principal Value
(D12) 0
although, I must confess, most probably it knows nothing about
finite-part integrals.
Dear Barton,
I feel you have already done an excellent work with the specfun
package and the additional packages you have prepared for Maxima,
which now tend to become standard packages and you should be
congratulated for your effort. I took the liberty above to suggest
some extensions according to my own small experience with ortho-
gonal polynomials, which, I feel, generally will not be a very difficult
task for you in case you can find the necessary time to proceed.
Please, forgive me for any incorrect comments above, which are
simply due to my being a mechanical-eletrical engineer without the
sufficient mathematical and computer algebra background. (Analogous
is the situation with Maxima as far as I am concerned.)
Many sincere thanks again for your precious extensions of Maxima.
With my best regards from Patras,
Nikos