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 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.