compare.jl

function quad42(f,a,b,tol; α=0.9, Θ=3.0)
    x, val, h = a, 0, (b-a)/10
    xx, f₀ = [x], f(x)
    err, evals = 0.0, 1
    
    while x < b
        fₘ, f₁ = f(x+h/2), f(x+h) # effektive Stufenzahl 2 (2 f-Auswertungen pro Schritt)
        evals += 2
        
        Q2 = h*(f₀+f₁)/2      # Trapezregel (Ordnung 2)
        Q4 = h*(f₀+4*fₘ+f₁)/6 # Simpson-Regel (Ordnung 4)
        
        if (ϵ = abs(Q4-Q2)) < tol
            err += ϵ
            x += h
            push!(xx,x)  # füge Knoten zum adaptiven Gitter hinzu
            f₀ = f₁      # FSAL-Trick
            val += Q4
        end
        
        h = min(0.9h*(tol/ϵ)^(1/3),3h,b-x) # adaptive Schrittweitenformel
    end
    return val, err, evals
end

function quad42₃(f,a,b,tol; α=0.9, Θ=3.0)
    x, val, h = a, 0, (b-a)/10
    xx, f₀ = [x], f(x)
    err, evals = 0.0, 1
    
    while x < b
        fₘ, f₁ = f(x+h/2), f(x+h) # effektive Stufenzahl 2 (2 f-Auswertungen pro Schritt)
        evals += 2
        
        Q2 = h*(f₀+f₁)/2      # Trapezregel (Ordnung 2)
        Q4 = h*(f₀+4*fₘ+f₁)/6 # Simpson-Regel (Ordnung 4)
        
        ϵₘ = abs(Q4-Q2)
        
        if (ϵ = ϵₘ ^ 2 / abs(Q4)) < tol
            err += ϵ
            x += h
            push!(xx,x)  # füge Knoten zum adaptiven Gitter hinzu
            f₀ = f₁      # FSAL-Trick
            val += Q4
        end
        
        h = min(0.9h*(tol/ϵ)^(1/3),3h,b-x) # adaptive Schrittweitenformel
    end
    return val, err, evals
end

using FastGaussQuadrature

function gauss_legendre(s)
    # für Interval [0,1]; Methode von Golub/Welsh (1969)
    k = 1:s-1
    β = k./sqrt((2*k-1).*(2*k+1));
    T = SymTridiagonal(zeros(s),β)
    λ, Q = eig(T)
    c = (1+λ)/2       # Knoten
    w = Q[1,:].^2     # Gewichte
    P = Q./(Q[1,:]')  # Matrix für interpolatorische Formeln
    return c,w,P
end

function embedded_gauss_legendre(m)   
    # eingebette Formeln der Stufenzahlen (4m+3,4m+2,2m) und Ordnungen (8m+3,4m+2,2m)
    
    global c, w, dw1, dw2
    s, s1, s2 = 4m+3, 4m+2, 2m
    w1, w2, e1 = zeros(s), zeros(s), zeros(s)
    e1[1] = 1
    
    ind1 = collect(1:s) 
    ind2 = collect(1:s)
    deleteat!(ind2,1:2:s)
    deleteat!(ind2,m+1)
    deleteat!(ind1,2m+2)
    
    c, w, P = gauss_legendre(s)
    k1 = length(ind1); w1[ind1] = P[1:k1,ind1]\e1[1:k1]
    k2 = length(ind2); w2[ind2] = P[1:k2,ind2]\e1[1:k2]
    
    dw1 = w-w1
    dw2 = w-w2
    
    return 2s
end

function quadrature_step(f,x,h)
    global nfcn
    fval = f.(x+h*c)
    evals = length(x)
    nfcn += length(c)
    err1 = h*dot(dw1,fval)
    err2 = h*dot(dw2,fval)
    return h*dot(w,fval), abs(err1)*(err1/err2)^2, evals # Wert der Quadraturformel, lokale Fehlerschätzung
end

function AdaptiveGaussQuadrature(f,a,b,tol,m)
    
    global nfcn = 0
    p = embedded_gauss_legendre(m)  
    
    x, val, err, h = a, 0, 0, (b-a)/10
    evals = 0
    xx = [x]
    while x < b
        Qp, ϵ, evs = quadrature_step(f,x,h)
        evals += evs
        if ϵ <= tol
            x += h
            push!(xx,x)  # füge Knoten zum adaptiven Gitter hinzu
            val += Qp
            err += ϵ
        end
        h = min(0.75h*(tol/ϵ)^(1/(p+1)),2h,b-x) # adaptive Schrittweitenformel
    end
    return val, err, evals
end

fns = [
    (-1.0,   1.0, (x) -> 1.0/(1e-4 + x^2)),
    ( 0.0,   1.0, (x) -> sqrt(x)*log(x)),
    (10.0, 110.0, (x) -> 2 + sin(3*cos(0.002 * (x-40)^2)))
]

Qs = [
    ("quad42", quad42),
    ("quad42₂", quad42₃),
    ("gauss₁", (f,a,b,tol) -> AdaptiveGaussQuadrature(f,a,b,tol,1)),
    ("gauss₂", (f,a,b,tol) -> AdaptiveGaussQuadrature(f,a,b,tol,2)),
    ("gauss₃", (f,a,b,tol) -> AdaptiveGaussQuadrature(f,a,b,tol,3))
]

function plot_it(fn, Qs)
    a,b,f = fn
    
    pows = 10.^(-1.0 .* (2:15))
    tols = reshape([pows 0.5.*pows]', 28)
    
    xscale("log")
    yscale("log")
    
    for Q in Qs
        name, I = Q
        vals = ((tol) ->
            let
                _,err,evals = I(f, a, b, tol)
                return 1.0*evals, err
            end).(tols)        
        
        vals = [p[i] for p in vals for i=1:2]
        vals = reshape(vals, (2, trunc(Int, length(vals)/2)))
            
        loglog(vals[1,:], vals[2,:], ".", label=name)
    end
    
    grid("on")
    xlabel("# f Auswertungen")
    ylabel("Fehler")
    legend()
end

plot_it(fns[1], Qs)
plot_it(fns[2], Qs)
plot_it(fns[3], Qs)