Hello Waldek, Thanks a lot for your response. Apologies for the delay in my reply. I will replace *Float *and not *DoubleFloat* and see if I'm still getting working solutions.
I will take a look into modifying my current code to modify it. As Tobias suggested, right now I also perhaps need to take a deep dive into various methods of numerical integration to get a hang of which one might be the best for my functions. best regards, Thejasvi On Thursday, April 15, 2021 at 12:05:29 AM UTC+2 Waldek Hebisch wrote: > 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/54ab9840-72ee-4263-bd0e-332c012cdaefn%40googlegroups.com.
