Alasdair McAndrew wrote:
>
> I have no particular attachment to Gauss-Kronrod integration - indeed
> Clenshaw-Curtis integration may well be easier to implement. My interest
> (as well as having a robust numeric integration routine in FriCAS) is
> seeing if an external numeric library can be used through FriCAS. And
> there are lots of these: QUADPACK, FFTW, ODE, ODEPACK to name but four
> Netlib ones. I don't see why any of us should have to reinvent the wheel
> when there are perfectly good open-source numeric routines already
> available. All I want is a FriCAS front end to them - and with arbitrary
> precision.
Well, we can use routines from other projects if they are
available. However, note that popular libraries typically
are limited to double precision. To get arbitrary precision
in most cases we will have to roll our own. The routine
I wrote is limited to double precision, but once we work
out algorithms we should be able to produce version for
arbitarary precision (as you note below we can produce
nodes and weights for Gaussian quadrature of arbitrary
order). In fact, from my point of view main reason
for using foreign libraries is speed -- assuming that
authors put more effort into speed optimizations that
we want to. However quadratures (and to some degree
ODE solvers) are special -- their execution time
should be dominated by user provided routines. Assuming
that we use good algorithm we should be at least as
fast as external libraries. Namely with good algorithm
we will have the number of calls to user routines and
calling FriCAS routine (I assume that user routine is
written in FriCAS) from FriCAS is faster than calling
it from outside of FriCAS.
In case of elliptic functions I do not know of any
open source system having routines as good as ours.
OTOH Lapack and Blas are definitely worth reusing
(we still should roll our owen arbitrary precision
versions but for machine precision Lapeck and Blas
have definite speed advantage).
So part important part of work is to identify good
routines -- I hope on some first hand recommendations
for good routines.
>
> For example, the nodes and weights above could be obtained with:
>
> digits(40)
> outputGeneral(40)
> x1 := solve(legendreP(9,x),1.0e-40)
> xs := [rhs(s) for s in x1]
> ws := [2/(1-s^2)/subst(D(legendreP(9,x),x),x=s)^2 for s in xs]
>
> And then you have to compute the Kronrod nodes, which interpolate between
> the Gaussian nodes...
--
Waldek Hebisch
--
You received this message because you are subscribed to the Google Groups
"FriCAS - computer algebra system" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To post to this group, send email to [email protected].
Visit this group at http://groups.google.com/group/fricas-devel.
For more options, visit https://groups.google.com/d/optout.
