Using GSL.jl integration routines in julia: integration_qawc

user_1117766 0
Viewed 222

Analogous to the example given in GSL.jl/examples/Quadrature.jl I am trying to integrate a function. However, since this function has a singularity, I need to use the cauchy weight. My idea was to use the following code

using GSL
function Q(p)
    ws_size = 200 
    ws     =  GSL.integration_workspace_alloc(ws_size)
    f_ = x -> 1/(x+p)
    f = GSL.@gsl_function(f_)
    result = Cdouble[0][1]
    epsrel = 1e-10 
    epsabs = 1e-10
    abserr = Cdouble[0][1]
    limit  = Csize_t[0][1]
    result = integration_qawc(f, 0., 1.e4, p, epsabs,epsrel,limit,ws,result,abserr)
    GSL.integration_workspace_free(ws)    
    return result
end

However, I get the following error

    UndefVarError: f_ not defined

    Stacktrace:
     [1] (::getfield(Main, Symbol("##117#118")))(::Float64, ::Ptr{Nothing}) at /home/varantir/.julia/packages/GSL/IVE5m   /src/manual_wrappers.jl:45
     [2] integration_qawc at /home/varantir/.julia/packages/GSL/IVE5m/src/gen/direct_wrappers/gsl_integration_h.jl:570 [inlined]
     [3] Q(::Float64) at ./In[250]:14

[4] top-level scope at In[251]:1

Which seems a little bit strange to me, since I clearly have defined f_. Any ideas?

gsl julia
1 Answers

for a weird reason, this doesn't throw errors but throws 0:

function Q(p)
    ws_size = 200 
    ws     =  GSL.integration_workspace_alloc(ws_size)
    f0(x::Float64)::Float64 = 1/(x+p)
    f = GSL.@gsl_function(f0)
    result = Cdouble[0][1]
    epsrel = 1e-10 
    epsabs = 1e-10
    abserr = Cdouble[0][1]
    limit  = Csize_t[0][1]
    result = integration_qawc(f, 0., 1.e4, p, epsabs,epsrel,limit,ws,result,abserr)
    GSL.integration_workspace_free(ws)    
    return result
end

From the docs of integration_qawc:

The adaptive bisection algorithm of QAG is used, with modifications to ensure that subdivisions do not occur at the singular point x = c. When a subinterval contains the point x = c or is close to it then a special 25-point modified Clenshaw-Curtis rule is used to control the singularity. Further away from the singularity the algorithm uses an ordinary 15-point Gauss-Kronrod integration rule.

Using an alternative, using QuadGK.jl:

using QuadGK
function G2(p)
    f(x)=1/(x+p)
    a = 0.0
    b = 1e4
    if a<-p<b
        res, err =  quadgk(f,a,-p,b,rtol=1e-10,atol=1e-10)
        return res
    else
        res, err =  quadgk(f,a,b,rtol=1e-10,atol=1e-10)
        return res
    end
end

from the QuadGK docs:

The algorithm is an adaptive Gauss-Kronrod integration technique: the integral in each interval is estimated using a Kronrod rule (2*order+1 points) and the error is estimated using an embedded Gauss rule (order points). The interval with the largest error is then subdivided into two intervals and the process is repeated until the desired error tolerance is achieved.

These quadrature rules work best for smooth functions within each interval, so if your function has a known discontinuity or other singularity, it is best to subdivide your interval to put the singularity at an endpoint. For example, if f has a discontinuity at x=0.7 and you want to integrate from 0 to 1, you should use quadgk(f, 0,0.7,1) to subdivide the interval at the point of discontinuity. The integrand is never evaluated exactly at the endpoints of the intervals, so it is possible to integrate functions that diverge at the endpoints as long as the singularity is integrable (for example, a log(x) or 1/sqrt(x) singularity).

The default order is 7, so is equivalent to the GSL integration.

user_121086990
Related