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.
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... On Sun, Oct 25, 2015 at 2:00 PM, Waldek Hebisch <hebi...@math.uni.wroc.pl> wrote: > Alasdair McAndrew wrote: > > > > The more I read about Lisp the more difficult it seems. It would > probably > > be just as easy for me to write a integrate.input file from scratch > > containing the Gauss-Kronrod code, > > Any reason for specifically wanting Gauss-Kronrod code? On > my todo list was writing adaptive Gauss routine. So I did > one using 9 point rule. Hopefuly it should be good enough > to get 14 digit accuracy with reasonable number of evaluation > points for large class of functions. This is preliminary > version, improvements welcome. > > BTW: I believe that adaptive quadratures, recursing on > subintervals are much better than nonadaptive. This one > is quite simplistic, in particular it can handle singularity > in sqrt, but fails on 1/sqrt. > > -------------<cut here>----------------- > > )abbrev package GAUSSI GaussianQuadrature > GaussianQuadrature : Exports == Implementation where > FT ==> DoubleFloat > Exports ==> with > adaptive : (FT -> FT, FT, FT, FT, > Integer, (FT -> FT, FT, FT) -> FT) -> FT > ++ adaptive(f, a, b, eps, n, g) is general adaptive > ++ quadrature on interval (a, b) using g as base > ++ (nonadaptive) rule and stopping when error is less than > ++ eps. adaptive signals error when number of recursion > ++ levels exceeds n. > ++ Note: agruments to g are f, midpoint of interval and > ++ half of length > gauss9 : (FT -> FT, FT, FT) -> FT > ++ gauss9(f, mid, h) computes 0 point Gaussian quadrature > ++ of f on the interval with midpoint mid and lenght 2*h > Implementation ==> add > > x1 := 0.3242534234_0380892904::DoubleFloat > x2 := 0.6133714327_0059039731::DoubleFloat > x3 := 0.8360311073_266357943::DoubleFloat > x4 := 0.9681602395_0762608983::DoubleFloat > w0 := 0.3302393550_0125976316::DoubleFloat > w1 := 0.3123470770_4000284007::DoubleFloat > w2 := 0.2606106964_0293546232::DoubleFloat > w3 := 0.1806481606_9485740405_8::DoubleFloat > w4 := 0.0812743883_6157441197_2::DoubleFloat > half := 0.5::DoubleFloat > > gauss9(f : FT -> FT, mid : FT, h : FT) : FT == > h1 := h*x1 > h2 := h*x2 > res := w0*f(mid) + w1*(f(mid + h1) + f(mid - h1)) > + w2*(f(mid + h2) + f(mid - h2)) > h3 := h*x3 > h4 := h*x4 > res := res + w3*(f(mid + h3) + f(mid - h3)) > + w4*(f(mid + h4) + f(mid - h4)) > h*res > > adaptive1(f : FT -> FT, mid : FT, h : FT, eps : FT, > n : Integer, g : (FT -> FT, FT, FT) -> FT, res0 : FT) : FT == > n < 0 => error "too many recursion levels" > h1 := half*h > mid1 := mid - h1 > mid2 := mid + h1 > res1 := g(f, mid1, h1) > res2 := g(f, mid2, h1) > nres := res1 + res2 > abs(res0 - nres) < eps => nres > eps1 := half*eps > adaptive1(f, mid1, h1, eps1, n - 1, g, res1) > + adaptive1(f, mid2, h1, eps1, n - 1, g, res2) > > adaptive(f : FT -> FT, a : FT, b : FT, eps : FT, > n : Integer, g : (FT -> FT, FT, FT) -> FT) : FT == > mid := half*(a + b) > h := half*(b - a) > res0 := g(f, mid, h) > adaptive1(f, mid, h, eps, n, g, res0) > > --------------------<cut here>------------------- > > -- > Waldek Hebisch > > -- > You received this message because you are subscribed to a topic in the > Google Groups "FriCAS - computer algebra system" group. > To unsubscribe from this topic, visit > https://groups.google.com/d/topic/fricas-devel/q18Av7P3jnM/unsubscribe. > To unsubscribe from this group and all its topics, send an email to > fricas-devel+unsubscr...@googlegroups.com. > To post to this group, send email to fricas-devel@googlegroups.com. > Visit this group at http://groups.google.com/group/fricas-devel. > For more options, visit https://groups.google.com/d/optout. > -- [image: http://www.facebook.com/alasdair.mcandrew] <http://www.facebook.com/alasdair.mcandrew> [image: https://plus.google.com/+AlasdairMcAndrew/posts] <https://plus.google.com/+AlasdairMcAndrew/posts> [image: https://www.linkedin.com/pub/alasdair-mcandrew/a/178/108] <https://www.linkedin.com/pub/alasdair-mcandrew/a/178/108> [image: https://twitter.com/amca01] <https://twitter.com/amca01> [image: http://numbersandshapes.net] <http://numbersandshapes.net> -- 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 fricas-devel+unsubscr...@googlegroups.com. To post to this group, send email to fricas-devel@googlegroups.com. Visit this group at http://groups.google.com/group/fricas-devel. For more options, visit https://groups.google.com/d/optout.
