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