On Wed, Apr 14, 2021 at 02:50:49PM +0200, thejasvi wrote:
> Hello Waldek and Ric
> Thanks so much for the function and comments you sent.
> The integration works pretty fast for the example function you showed.
>
> For my own use-case, I expect to need very high precision (eg. the
> Mathematica code I'm trying to port uses 300 digits for the same model).
> Is DoubleFloat the correct Type to use?
No, DoubleFloat is limited to about 15 digits. You need to use Float
and use something like 'digits(300)' to set precision. Note: if you
really need 300 digits of final precision, you need to set bigger
precision for intermediate calculation.
>
> Given that the target is 300 digits precision, I tried to run the adaptive
> integration with lower and lower errors (< 1.0e-16).
> It ended up taking a much longer time (didn't complete in 3mins). Do any of
> you have ideas on how to overcome this.
This routine is for DoubleFloat accuracy, using it as-is you can not
get better than what DoubleFloat allows. For better accuracy you
need the following:
- replace DoubleFloat by Float.
- provide pl and wl of sufficient accuracy (that is 300 digits for
300 results)
- probably increase order: current version uses Gaussian kwadrature
on 9 point, which gives order 17. At first glance 75 points
looks like reasonable choice for 300 digits accuracy.
pl are nodes of Gaussian quadrature, they are zeros of appropriate
Legendre polynomial. wl are weights. My code assumes that number
of points is odd, then 0 is one of nodes, other nodes are symmetric.
I took adantage of symmetry to remember only half of nodes and
woights. Below is code that I used to compute nodes and weights
(this one for 45 points):
digits(70)
pn := legendreP(45, x)
sols := solve(%, 1.0e-60)
solp := map(x +-> retract(rhs(x))@Float, last(sols, 23))
dpn := D(pn, x)
wpp := [2/((1 - xi^2)*eval(dpn, x = xi)^2) for xi in solp]
wlf := map(retract, wpp)@List(Float)
plf := rest(solp)
wl := wlf::List(DoubleFloat)
pl := plf::List(DoubleFloat)
Of course, for Float the last two lines should be replaced by
assignments. Note that 'solve' needs substantial time: it
uses robust method which however is slower than Newton when
you need large accuracy. So, using 'solve' with limited
accuracy and then Newton to improve approximation may be
faster. Extra remark: you should compute nodes and weights
once and then use for all integrations with given accuracy.
Extra remarks:
- For highly oscilatory integrals you should use specialized
method. I did not look carefuly enough at your code to know
if you can use general method (Gaussian quadrature) or
need special one
- ATM in FriCAS Bessel functions have only DoubleFloat accuracy.
Bessel functions of half-integer order are elementary, so
for fixed (and hopefuly moderate) order you can use explicit
formulas to provide your own version
--
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 view this discussion on the web visit
https://groups.google.com/d/msgid/fricas-devel/20210414220521.GA16203%40math.uni.wroc.pl.
