juliamath / quadgk.jl Goto Github PK
View Code? Open in Web Editor NEWadaptive 1d numerical Gauss–Kronrod integration in Julia
License: MIT License
adaptive 1d numerical Gauss–Kronrod integration in Julia
License: MIT License
Suggestion to make it more clear in the documentation that a function is required in QuadGK to evaluate the data and that
for purely numerical input (e.g., as recorded by instruments), other packages are recommended (e.g. Romberg.jl for equally spaced data or Trapz.jl for arbitrarily spaced data).
This issue is used to trigger TagBot; feel free to unsubscribe.
If you haven't already, you should update your TagBot.yml
to include issue comment triggers.
Please see this post on Discourse for instructions and more details.
If you'd like for me to do this for you, comment TagBot fix
on this issue.
I'll open a PR within a few hours, please be patient!
While this package allows evaluation of integrands which return arrays, it's not efficient for the integrand to be an allocating f(t)
. Instead, quadgk
should support a non-allocating form f!(out,t)
for this case, for this version of the function do minimal internal allocation by caching necessary arrays and doing everything in-place.
Hi there, I would like to ask a question regarding quadrature points and weights in the intervals [0, +Inf]
and [-Inf, +Inf]
. At the moment, it seems that I could not get quadrature points and weights in these intervals using the command gauss(N,a,b)
. However, it is possible to find quadgk(x -> exp(-x^2), -Inf, +Inf)
or quadgk(x -> exp(-x^2), 0, +Inf)
. So, my question here is that will the command gauss(N,a,b)
support the intervals [0, Inf]
and [-Inf, +Inf]
in a near future? Or I have to perform a change of variable to scale these intervals into finite intervals to get the quadrature points and weights?
Thank you for a great package!
I would like to use Zygote on quadgk
but below is a self-contained example, where it does not work because quadgk
calls setindex!
. Can I do anything to work around this? I coded up a small example of differentiating quadgk
of a constant function and it did work so it seems that there is hope for Zygote and QuadGK to play together nicely.
using Zygote, QuadGK
"P(0<X<x)"
function F(c, x)
if x < c
return 0.
elseif x < 1.
return exp(c)*(x-c)^2/2
elseif x < 1. + c
return exp(c)*(1-c)*(2x - c - 1.)/2.
elseif x < 2.
return exp(c)*(4x - x^2 - 2c - 2)/2.
else
return exp(c)*(1-c)
end
end
"pdf of X"
function f(c, x)
if x < c || x > 2.
return 0.
elseif x < 1.
return exp(c)*(x-c)
elseif x < 1. + c
return exp(c)*(1-c)
else
return exp(c)*(2-x)
end
end
function D(c)
exp(c)*(1-c)
end
function W(c1, c2)
val, acc = quadgk(t -> F(c1, t)*f(c2, t), c1, 2)
val
end
function prob_of_win(c1, c2)
((1-D(c1))*D(c2) + W(c2, c1))/(D(c1) + D(c2) - D(c1)*D(c2))
end
gradient(prob_of_win, 0.25, 0.25) # ERROR: Mutating arrays is not supported -- called setindex!(Vector{QuadGK.Segment{Float64, Float64, Float64}}, ...)
rulekey(::Type{BigFloat}, n) = (BigFloat, precision(BigFloat), n)
rulekey(T,n) = (T,n)
The function returns either a triplet or a pair, this would make it easier to reason about if it returned always the same structure and types
@stevengj Do you have any input on what should be done here before this is tagged and registered?
Following is my code to solve two-dimensional integral in Julia.
using Statistics,Distributions,Distributed,StatsFuns,Roots,QuadGK
p=3;n0=3;n1=1;n2=4;d=0.25;
inARL=200;n=n1+n2;w1=chisqinvcdf(p,1-(n0-n1)/n2);inpa1=chisqcdf(p,w1);
function inpa(h2)
function fvv(v)
return(nchisqcdf(p,(n1/n2)v,(n*h2)/n2)*chisqpdf(p,v));end;
return(1/(1-inpa1-quadgk(fvv,w1,Inf,rtol=1e-10)[1])-inARL);end;
h=fzero(inpa,0.25,50);
outpa1=nchisqcdf(p,n1*d^2,w1);
function fz1(z)
function fv1(v)
lam=(n1/n2)*((z+sqrt(n1)*d)^2+v)+(2*n*sqrt(n1)*d*z)/(n2)+((n*d)^2)/(n2);
return(nchisqcdf(p,lam,(n*h)/n2)*chisqpdf(p-1,v));end;
return(normpdf(z)*quadgk(fv1,0,Inf,rtol=1e-10)[1]);end;
int1=quadgk(fz1,-Inf,-sqrt(w1)-sqrt(n1)*d,rtol=1e-10)[1];# problem
An issue was caught in SciML/DataInterpolations.jl#178 where tests started failing with [email protected]
I have made a MWE to illustrate the issue (using a simple LinearInterpolation from DataInterpolations):
u = collect(0:10)
t = collect(0:10)
A = LinearInterpolation(u,t)
qint, err = quadgk(A, 0.0, 10.0)
It is a straight line with slope 1. So the integral should be 50.0.
with [email protected]
julia> qint, err = quadgk(A, 0.0, 10.0)
(50.0, 0.0)
with [email protected]
julia> qint, err = quadgk(A, 0.0, 10.0)
(51.76911234248546, 4.274582013885939)
I just started experimenting with the recent DynamicQuantities.jl package, which aims to provide the same functionality as Unitful.jl but with a statically-typed Quantity
to improve performance. It looks like QuadGK
currently doesn't support these DynamicQuantities Quantity
types, but I'm not familiar enough with the codebase to quite understand why.
I'm interested in contributing by generating a PR for this, but would really appreciate a nudge in the right direction. It looks like the issue spawns from approximately here, where either a @generated
cachedrule
isn't being produced at compile-time or maybe the generic cachedrule
below it isn't working.
Any ideas that might help me understand what I can do to help integrate this functionality?
Thanks,
Mike
Minimal reproduction (Julia v1.9.3):
julia> using DynamicQuantities # [06fc5a27] DynamicQuantities v0.7.0
julia> using QuadGK # [1fd47b50] QuadGK v2.8.2
julia> quadgk(t -> 5u"m/s"*t, 0u"s", 10u"s")
ERROR: MethodError: no method matching cachedrule(::Type{Quantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}}}, ::Int64)
Closest candidates are:
cachedrule(::Union{Type{Complex{BigFloat}}, Type{BigFloat}}, ::Integer)
@ QuadGK C:\Users\*\.julia\packages\QuadGK\XqIlh\src\gausskronrod.jl:253
cachedrule(::Type{T}, ::Integer) where T<:Number
@ QuadGK C:\Users\*\.julia\packages\QuadGK\XqIlh\src\gausskronrod.jl:264
quadgk
is not compatible with AD (well, with ForwardDiff.jl) because of cachedrule
. I am trying to solve an optimization problem using JuMP where some nonlinear constraints involve integrals. I can compute them using other methods, but it would be nice if QuadGK "just worked"
Here is the rather uninformative error
ERROR: StackOverflowError:
Stacktrace:
[1] cachedrule(#unused#::Type{ForwardDiff.Dual{ForwardDiff.Tag{JuMP.var"#137#139"{typeof(integratesegmentgk)}, Float64}, Float64, 5}}, n::Int64) (repeats 79984 times)
@ QuadGK C:\Users\beasont\.julia\packages\QuadGK\ENhXl\src\gausskronrod.jl:253
The reference to integratesegmentgk
is my user-defined function in the constraint. It simply calls quadgk
, but for completeness:
function integratesegmentgk(mL,mU,a,b,c)
hfun(x) = exp(-a -b*x -c*x^2)
res,tol = quadgk(hfun,mL,mU)
return res
end
I am using FastGaussQuadrature.jl in the meantime.
I used to be able to do
julia> g(y) = quadgk(x -> x, y, 1.0)
g (generic function with 1 method)
julia> g(0.0)
(0.5, 0.0)
julia> g(1.0)
(0.0, 0.0)
Doing an integral where the upper and lower bound were identical gave the answer 0.0, as expected. After upgrading to v2.9.2, I get the error:
julia> g(1.0)
ERROR: DomainError with 1.0:
roundoff error detected near endpoint of the initial interval (1.0, 1.0)
Stacktrace:
[1] throw_endpoint_error(x::Float64, a::Float64, b::Float64)
@ QuadGK ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:186
[2] #check_endpoint_roundoff#48
@ ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:181 [inlined]
[3] check_endpoint_roundoff
@ ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:179 [inlined]
[4] #6
@ ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:15 [inlined]
[5] ntuple
@ ./ntuple.jl:48 [inlined]
[6] do_quadgk(f::var"#5#6", s::Tuple{…}, n::Int64, atol::Nothing, rtol::Nothing, maxevals::Int64, nrm::typeof(LinearAlgebra.norm), segbuf::Nothing)
@ QuadGK ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:13
[7] #51
@ ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:257 [inlined]
[8] handle_infinities(workfunc::QuadGK.var"#51#52"{…}, f::var"#5#6", s::Tuple{…})
@ QuadGK ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:146
[9] #quadgk#50
@ QuadGK ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:256 [inlined]
[10] quadgk
@ QuadGK ~/.julia/packages/QuadGK/PWpum/src/adapt.jl:254 [inlined]
[11] g(y::Float64)
@ Main ./REPL[2]:1
[12] top-level scope
@ REPL[3]:1
which seems to be caused by #91 .
Is there anyway to restore the functionality to accept identical lower and upper bounds? The above function is just an MWE, but it's not uncommon to have an integral inside a function where the lower bound is an argument and can vary to the upper bound.
DataStructures 0.18 is out, can this package be revised to support it? I am having conflicts between this package and the latest version of DiffEqBase due to this dependency version.
When I run the integral of f(x)="0.21702437*(0.6227007 + 1.4335294/(x + 1.6759317))/(x + 1.52924537 - 0.2316349/x)" from x=0 to x=1, by running this code in Python:
jl.eval("""
my_tuple = quadgk(x -> 0.21702437*(0.6227007 + 1.4335294/(x + 1.6759317))/(x + 1.52924537 - 0.2316349/x), 0, 1, rtol=0.01)
println(my_tuple)
""")
I consistently get a value of "0.17259607425281273" for my_tuple[1].
I would expect to get a value of infinity or undefined.
Interestingly, when I use PySR python package's eval_tree_array and I pass in a tree whose definition is:
x1 = Node(;feature=1)
tree = 0.21702437*(0.6227007 + 1.4335294/(x1 + 1.6759317))/(x1 + 1.52924537 - 0.2316349/x1)
tree = Node{Float32}(tree)
my_tuple = quadgk(x -> (eval_tree_array(tree, reshape([Float32(x)], 1, 1), options))[1][1], 0, 1, rtol=0.01)
I consistently get a value of "0.14588867513106551" for my_tuple[1].
Perhaps the difference between 0.17259 and 0.14588 is because I converted the tree to Float32 (instead of Float64).
If you open up Desmos.com and plot this function f(x) over x=0 to x=1, you will find that the function goes to negative infinity as x approaches 0.16 from the to-the-right direction, and the function goes to positive nfinity as x approaches 0.16 from the to-the-left direction.
And I'm pretty sure that in reality, the integral of this function is undefined. Although perhaps, the negative infinity cancels out the positive infinity, which means the value of the definite integral should be finite?
Thanks for all the help, in advance.
I ran
@code_warntype quadgk( x->x, 0., 1. )
today and it says, that the return type doesn't get inferred. Is it me or is it a bug?
The issue is that segbufs don't always work with infinite limits, such as when the limits have units.
This is a MWE
using QuadGK, Unitful
Tu = typeof(1.0u"m")
segbuf = QuadGK.alloc_segbuf(Tu, Tu, Tu)
quadgk(x -> 1/(1+oneunit(x)/x), 0.0u"m", 1.0u"m", segbuf=segbuf) # OK (0.3068528194400547 m, 1.9931833961095435e-11 m)
quadgk(x -> 1/(1+oneunit(x)/x), 0.0u"m", Inf*u"m", segbuf=segbuf) # broken
which gives the error
ERROR: DimensionError: m and 0.0 are not dimensionally compatible.
Stacktrace:
[1] convert(#unused#::Type{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}, x::Float64)
@ Unitful ~/.julia/packages/Unitful/jiUeO/src/conversion.jl:106
[2] QuadGK.Segment{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}(a::Float64, b::Float64, I::Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, E::Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}})
@ QuadGK ~/.julia/dev/QuadGK/src/evalrule.jl:3
[3] convert(#unused#::Type{QuadGK.Segment{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}, s::QuadGK.Segment{Float64, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}})
@ QuadGK ~/.julia/dev/QuadGK/src/evalrule.jl:10
[4] setindex!(A::Vector{QuadGK.Segment{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}, x::QuadGK.Segment{Float64, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}, i1::Int64)
@ Base ./array.jl:969
[5] macro expansion
@ ./broadcast.jl:974 [inlined]
[6] macro expansion
@ ./simdloop.jl:77 [inlined]
[7] copyto!
@ ./broadcast.jl:973 [inlined]
[8] copyto!
@ ./broadcast.jl:926 [inlined]
[9] materialize!
@ ./broadcast.jl:884 [inlined]
[10] materialize!
@ ./broadcast.jl:881 [inlined]
[11] do_quadgk(f::QuadGK.var"#18#27"{var"#11#12", Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}, s::Tuple{Float64, Float64}, n::Int64, atol::Nothing, rtol::Nothing, maxevals::Int64, nrm::typeof(LinearAlgebra.norm), segbuf::Vector{QuadGK.Segment{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}})
@ QuadGK ~/.julia/dev/QuadGK/src/adapt.jl:39
[12] #50
@ ~/.julia/dev/QuadGK/src/adapt.jl:235 [inlined]
[13] handle_infinities(workfunc::QuadGK.var"#50#51"{Nothing, Nothing, Int64, Int64, typeof(LinearAlgebra.norm), Vector{QuadGK.Segment{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}}, f::var"#11#12", s::Tuple{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}})
@ QuadGK ~/.julia/dev/QuadGK/src/adapt.jl:126
[14] quadgk(::Function, ::Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, ::Vararg{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}; atol::Nothing, rtol::Nothing, maxevals::Int64, order::Int64, norm::Function, segbuf::Vector{QuadGK.Segment{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}})
@ QuadGK ~/.julia/dev/QuadGK/src/adapt.jl:234
[15] top-level scope
@ REPL[11]:1
This is possibly related to #30, but not sure. I am integrating a function that is basically the Gamma PDF (very spiked near zero, but tapers off as x -> infinity). The function does not return NaN for any value (I fix gpdf(0)=0, or else it would give NaN here).
quadgk(x->gpdf(x,l1,d1,s1),0,Inf)
ERROR: DomainError with 0.9999999999999964:
integrand produced NaN in the interval (0.9999999999999929, 1.0)
Stacktrace:
[1] evalrule(::getfield(QuadGK, Symbol("##12#18")){getfield(Main, Symbol("##41#42")),Float64}, ::Float64, ::Float64, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::typeof(LinearAlgebra.norm)) at C:\Users\tbeason\.julia\packages\QuadGK\v6Zrs\src\evalrule.jl:41
[2] adapt(::Function, ::Array{QuadGK.Segment{Float64,Float64,Float64},1}, ::Float64, ::Float64, ::Int64, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Int64, ::Float64, ::Float64, ::Int64, ::typeof(LinearAlgebra.norm)) at C:\Users\tbeason\.julia\packages\QuadGK\v6Zrs\src\adapt.jl:39
[3] do_quadgk(::getfield(QuadGK, Symbol("##12#18")){getfield(Main, Symbol("##41#42")),Float64}, ::Tuple{Float64,Float64}, ::Int64, ::Nothing, ::Nothing, ::Int64, ::typeof(LinearAlgebra.norm)) at C:\Users\tbeason\.julia\packages\QuadGK\v6Zrs\src\adapt.jl:28
[4] handle_infinities(::getfield(Main, Symbol("##41#42")), ::Tuple{Float64,Float64}, ::Int64, ::Nothing, ::Nothing, ::Int64, ::Function) at C:\Users\tbeason\.julia\packages\QuadGK\v6Zrs\src\adapt.jl:83
[5] #quadgk#21(::Nothing, ::Nothing, ::Int64, ::Int64, ::Function, ::typeof(quadgk), ::Function, ::Float64, ::Float64) at C:\Users\tbeason\.julia\packages\QuadGK\v6Zrs\src\adapt.jl:155
[6] #quadgk#20 at C:\Users\tbeason\.julia\packages\QuadGK\v6Zrs\src\adapt.jl:155 [inlined]
[7] quadgk(::Function, ::Int64, ::Float64) at C:\Users\tbeason\.julia\packages\QuadGK\v6Zrs\src\adapt.jl:151
[8] top-level scope at REPL[109]:1
The error indicates that the NaN arises in the handle_infinities
function, which I guess it where you transform the integration domain from (0,Inf) to something usable. If I use a large number instead of Inf, I do not get the error.
MWE
using SpecialFunctions, QuadGK
l1 = exp(-1.2604534253)
d1 = -0.916695708
s1 = 8.428294026
function gpdf(t,lambda,delta,sigma)
t == 0 && return 0.0
d2 = delta^(-2)
constantpiece = lambda*abs(delta/sigma)/gamma(d2)
LT = lambda*t
varpiece = LT^(1/(delta*sigma)-1) * exp(-LT^(delta/sigma))
val = constantpiece*varpiece
return val
end
quadgk(x->gpdf(x,l1,d1,s1),0,Inf)
quadgk(x->gpdf(x,l1,d1,s1),0,1e30)
Since a few days I notice below issue (both Julia 0.6 and Julia 0.7-):
ERROR: QuadGK: fetch failed to get commit 064d2f0, please file an issue at https://github.com/JuliaMath/QuadGK.jl/issues
Stacktrace:
[1] fetch(::Base.LibGit2.GitRepo, ::String, ::String) at ./pkg/write.jl:25
[2] #3 at ./pkg/write.jl:53 [inlined]
[3] with(::Base.Pkg.Write.##3#4{String,String}, ::Base.LibGit2.GitRepo) at ./libgit2/types.jl:608
[4] resolve(::Dict{String,Base.Pkg.Types.VersionSet}, ::Dict{String,Dict{VersionNumber,Base.Pkg.Types.Available}}, ::Dict{String,Tuple{VersionNumber,Bool}}, ::Dict{String,Base.Pkg.Types.Fixed}, ::Dict{String,VersionNumber}, ::Set{String}) at ./pkg/entry.jl:542
[5] update(::String, ::Set{String}) at ./pkg/entry.jl:461
[6] (::Base.Pkg.Dir.##4#7{Array{Any,1},Base.Pkg.Entry.#update,Tuple{String,Set{String}}})() at ./pkg/dir.jl:36
[7] cd(::Base.Pkg.Dir.##4#7{Array{Any,1},Base.Pkg.Entry.#update,Tuple{String,Set{String}}}, ::String) at ./file.jl:70
[8] #cd#1(::Array{Any,1}, ::Function, ::Function, ::String, ::Vararg{Any,N} where N) at ./pkg/dir.jl:36
[9] update() at ./pkg/pkg.jl:228
As mentioned on discourse, the kronrod
routine yields NaN for larger n
, e.g.
kronrod(550,1e-11,40)
It looks like this is due to spurious underflow in the _kronrodjacobi
function. It forms arrays s
and t
by calculations that involve repeatedly multiplying by the b
matrix, corresponding to the squared off-diagonal elements of the Jacobi matrix. For the canonical interval (-1,1) with a constant weight function, the entries of the b
matrix asymptote to 0.25, so this effectively divides by 4 over and over and eventually underflows.
The underflow is spurious because the final Kronrod–Jacobi matrix only depends on ratios of s
and t
entries, so we are free to rescale them by any overall constant factor that we want. This fact is actually mentioned in a somewhat obscure way by Laurie's paper:
Actually, as recommended by Reichel [23], the coefficients were scaled to the interval (−2, 2) before performing the computation, and the resulting points and weights (obtained by Gautschi’s routine gauss from [10]) were scaled back to (−1, 1). On a machine with a small exponent range, this subterfuge is necessary to avoid underflow of the modified moments.
Rescaling the interval to (-2,2) works because it multiplies the integrals by 2, and hence the squared Jacobi matrix entries get multiplied by 4, causing b
to asymptote to 1 and removing the exponential shrinkage.
We could do the same trick, but it's not clear how to generalize this to arbitrary Jacobi matrices. Instead, it should suffice to normalize the s
and t
matrices (e.g. scaling by their infinity norm) after each iteration. This should only slow things down by a small constant factor (1–2% in my current experiments).
Hi there,
I would like to ask you a question regarding quadrature rules for the weight function w(x)=1/(1+x^2)
. I hope this is fine to post my question here. My question is that I would like to compute the integral I = \int_{a}^{b} w(x)f(x)dx
using quadrature rules, where w(x)
is a weight function and f(x) can be approximated by a polynomial of degree 2N−1
or less on [a, b]
. This is because it is costly to compute f(x). Therefore, I am looking for a quadrature rule for the weight function w(x)=1/(1+x^2)
. I had a look at https://en.wikipedia.org/wiki/Gaussian_quadrature. However, I could not find a quadrature rule for the weight function that I need. But I found that your package can compute the integral I for arbitrary weight functions. So, do you think that it is fine that I use the function "QuadGK.gauss" to obtain the quadrature points x[i]
and weights w[i]
and can compute the integral I from sum(w .* f.(x))
. Or do you have any suggestions for me? Thank you very much in advance!
Right now the I and E types (returned from QuadGK) seem undefined in Segment
:
# integration segment (a,b), estimated integral I, and estimated error E
struct Segment
a::Number
b::Number
I
E
end
isless(i::Segment, j::Segment) = isless(i.E, j.E)
It could be replaced by something like:
struct Segment{T<:Real}
a::Number
b::Number
I::T
E::T
end
isless(i::Segment, j::Segment) = isless(i.E, j.E)
Two minimal working examples:
using QuadGK, Unitful
f(x) = x^2u"m"
quadgk(f, 0, 1)
and
using QuadGK, Unitful
f(x) = x^2
quadgk(f, 0u"m", 1u"m")
both crash Julia for me.
Seems to be directly related to #13. In hindsight, this probably could have been a comment in that thread. Should I remove this issue and move its contents there?
Hi,
Is there another way than using Sympy from Python?
If i do need to use it, as my integral is within a function and part of a larger numerical problem, should i then reconvert all my variables into Julia types (I mean versus symbolic - type Sym- that I had to convert them into to use integrate)?
Many thanks.
Could you please tag a new release for quadgk!?
Hi, I want to cite this package in a book and can't figure out the author field. So far I'm using
@misc{QuadGK,
title="{QuadGK.jl}",
year=2021,
note="Julia Package",
howpublished="https://github.com/JuliaMath/QuadGK.jl"
}
which works well with the publisher's LaTeX style. So .... author = ?
Thanks,
-- Tim
This seems wrong:
julia> quadgk_count(x -> 1, 0, 1, order=1)
(2.4444444444019076, 0.5555555555356156, 23333345)
QuadGK is holding my DataStructures to v0.15.
What is the reason for these lines in the Project.toml?
[compat]
DataStructures = "0.11, 0.12, 0.13, 0.14, 0.15"
Following is my code to solve two-dimensional integral in Julia.
using Statistics,Distributions,Distributed,StatsFuns,Roots,QuadGK
p=3;n0=3;n1=1;n2=29;d=0.25;
inARL=200;n=n1+n2;w1=chisqinvcdf(p,1-(n0-n1)/n2);inpa1=chisqcdf(p,w1);
function inpa(h2)
function fvv(v)
return(nchisqcdf(p,(n1/n2)v,(n*h2)/n2)*chisqpdf(p,v));end;
return(1/(1-inpa1-quadgk(fvv,w1,Inf,rtol=1e-10)[1])-inARL);end;
h=fzero(inpa,0.25,50);
outpa1=nchisqcdf(p,n1*d^2,w1);
function fz1(z)
function fv1(v)
lam=(n1/n2)*((z+sqrt(n1)*d)^2+v)+(2*n*sqrt(n1)*d*z)/(n2)+((n*d)^2)/(n2);
return(nchisqcdf(p,lam,(n*h)/n2)*chisqpdf(p-1,v));end;
return(normpdf(z)*quadgk(fv1,0,Inf,rtol=1e-10)[1]);end;
int1=quadgk(fz1,-Inf,-sqrt(w1)-sqrt(n1)*d,rtol=1e-10)[1];# problem
JuliaError: Exception 'UndefVarError: integral
not defined' occurred while calling julia code:
When I ran this code:
"
using QuadGK
try
integral, err = quadgk(x -> exp(1/(x-0.5)), 0, 1, rtol=1e-8)
catch
print("Error")
end
print(integral)
print(err)
"
I am generating a bunch of expression trees somewhat randomly, where each expression tree represents a particular function of x.
For example, one expression tree could be:
*
/
cos 5
|
x0
which represents the function "cos(x0)*5". I can call this function tree_as_func(x).
I need to run "quadgk(x -> tree_as_func(x), 0, 1, rtol=1e-8)" where "tree_as_func" is the function of x that the current expression tree represents.
I do not know if "tree_as_func(x)" blows up to infinity on the domain x=[0, 1) or if its integral over x=[0, 1) is infinity.
As a result, I need to be able to catch errors. But this try-catch block fails to catch the error.
How can I do this?
Hi,
wanted to know if this is the expected behavior. When doing
function f!(b, x)
b[1] = exp(-x^2)
b[2] = exp(-x^2)
return nothing
end
buf = zeros(Float64, 2); quadgk!(f!, buf, 0.0, Inf)
I get ERROR: MethodError: objects of type QuadGK.InplaceIntegrand{typeof(f!),Array{Float64,1},Array{Float64,1}} are not callable
. For finite upper bounds though everything seems to work. In my use case however, I need to compute integrals for which the upper bound is Inf, and I do not a priori know if and when I can truncate the integral for some large upper bound.
Dear QuadGK developers:
I'm attempting to calculate multiple 1D integrals in parallel using GPU via CUDA.jl. However currently QuadGK does not seem to be GPU-compatible (with CUDA.jl). A simple example below generates errors ( see the attached text file ). Is it possible, perhaps for the 1D integrals with limited options, to become compatible with CUDA.jl? Thanks.
using CUDA
using QuadGK
Z = Array([z for z in 1:1:4]);
Zcu = CuArray([z for z in 1:1:4]);
function f_int(z_ob)
iii(z) = z_ob - z
return quadgk(z -> iii(z), 0, 1, rtol=1e-4)[1]
end
f_int.(Z) # This works
f_int.(Zcu) # This does not work
When you follow the deprecation advice for using this package on Julia 0.6 and beyond,
julia> quadgk()
ERROR: quadgk() has been moved to the package QuadGK.jl.
Run Pkg.add("QuadGK") to install QuadGK on Julia v0.6 and later, and then run `using QuadGK`.
It does not actually work due to a clash with the deprecation warning itself.
julia> using QuadGK
julia> quadgk()
WARNING: both QuadGK and Base export "quadgk"; uses of it in module Main must be qualified
ERROR: UndefVarError: quadgk not defined
I believe the above was working on MASTER sometime over Winter 2016 / Spring 2017, as I had added using QuadGK
to all my codes to make the deprecation warning go away!
Instead you must import the function(s) you want explicitly,
julia> import QuadGK.quadgk
julia> quadgk
quadgk (generic function with 3 methods)
Or reference the QuadGK package every time you call the function,
julia> using QuadGK
julia> QuadGK.quadgk
quadgk (generic function with 3 methods)
See also discussion at:
JuliaLang/julia#19741
Can the QuadGK functions be forcibly exported into the namespace when the packaged is using
'd ? That would seem to be the expected behaviour by the user.
Hi, I have just tried to integrate a Laplace distribution as
using QuadGK, Distributions
pc = Laplace(0, 1)
quadgk(x -> pdf(pc, x), -1e8, 1e8, rtol=1e-8)
(0.0, 0.0)
I obtain zero, which is obviously wrong as it should be one. The manual says that I should add an explicit point, when the point is non-smooth, which I do
quadgk(x -> pdf(pc, x), -1e8, 0, 1e8, rtol=1e-8)
(0.0, 0.0)
But it does not help either. I think that the method just omit to sample around the zero sufficiently.
A simple syntax error when calling QuadGK returns an error giving the impression that QuadGK is not installed:
julia> using QuadGK
julia> g(x)=x^2
g (generic function with 1 method)
julia> quadgk(f,(0,1))
ERROR: quadgk(f, (0, 1)) has been moved to the package QuadGK.jl.
Run Pkg.add("QuadGK") to install QuadGK on Julia v0.6 and later, and then run `using QuadGK`.
Stacktrace:
[1] #quadgk#771(::Array{Any,1}, ::Function, ::Function, ::Vararg{Any,N} where N) at ./deprecated.jl:698
[2] quadgk(::Function, ::Tuple{Int64,Int64}) at ./deprecated.jl:698
julia>
This can cause quite a bit confusion for students learning Julia (and even some of the more familiar!).
By the way, great package!
I have a function if I have a function f(x)
which is the integral another function (using QuadGK) and x
is the upper limit of integration. If I want to take an automatic derivative of that function using ForwardDiff, it results in a stack overflow. Here are instructions to reproduce the error
using QuadGK
using ForwardDiff
D(f::Function, x::Number) = ForwardDiff.derivative(f, x)
f(x) = quadgk(y -> y^2, 0, x)
D(x -> quadgk(y -> y^2, 0, x), 1)
ERROR: StackOverflowError:
Stacktrace:
[1] #quadgk#15(::Array{Any,1}, ::Function, ::Function, ::ForwardDiff.Dual{ForwardDiff.Tag{##3#5,Int64},Float64,1}, ::ForwardDiff.Dual{ForwardDiff.Tag{##3#5,Int64},Float64,1}) at /Users/mason/.julia/v0.6/QuadGK/src/QuadGK.jl:241 (repeats 65479 times)
[2] derivative(::##3#5, ::Int64) at /Users/mason/.julia/v0.6/ForwardDiff/src/derivative.jl:14
[3] D(::Function, ::Int64) at ./REPL[4]:1
Some problem with the DOCUMENTER_KEY
ssh key according to the logs.
Initialized empty Git repository in /tmp/jl_B4CaCu/.git/
[email protected]: Permission denied (publickey).
fatal: Could not read from remote repository.
Please make sure you have the correct access rights
and the repository exists.
┌ Error: Git failed to fetch [email protected]:JuliaMath/QuadGK.jl.git
│ This can be caused by a DOCUMENTER_KEY variable that is not correctly set up.
│ Make sure that the environment variable is properly set up as a Base64-encoded string
│ of the SSH private key. You may need to re-generate the keys with DocumenterTools.
└ @ Documenter ~/.julia/packages/Documenter/6vUwN/src/Documenter.jl:552
I think there's something wrong with the way quadgk handles this integral, the error estimate is almost zero but the answer is substantially inaccurate. Just replacing 0.499 with 0.4999 fixes it:
julia> check(f,a,b) = let (x1,e1) = QuadGK.quadgk(f,a,b); (x2,e2) = QuadGK.quadgk(f,big(a),big(b)); (Float64(abs(x1-x2)), e1); end
check (generic function with 1 method)
julia> check(x -> exp(abs(x-0.499)), 0, 1)
(1.0000000833402914e-6, 2.220446049250313e-16)
julia> check(x -> exp(abs(x-0.4999)), 0, 1)
(7.284892445276705e-17, 1.1102230246251565e-16)
It's an easy integral, so zero error estimate can be expected, but there could be something going wrong with the integral value it actually returns.
This is probably related to #38, I am trying to integrate with quadgk a function from 1 to infinity which should be converging, and encounter this error (julia 1.8.5 and QuadGK v2.8.2).
using QuadGK
quadgk(x->1/x^1.1,1,Inf)
DomainError with 0.9999999999999964:
integrand produced NaN in the interval (0.9999999999999929, 1.0)
Stacktrace:
[1] evalrule(f::QuadGK.var"#14#23"{var"#1#2", Float64, Float64}, a::Float64, b::Float64, x::Vector{Float64}, w::Vector{Float64}, gw::Vector{Float64}, nrm::typeof(LinearAlgebra.norm))
@ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/evalrule.jl:37
[2] adapt(f::QuadGK.var"#14#23"{var"#1#2", Float64, Float64}, segs::Vector{QuadGK.Segment{Float64, Float64, Float64}}, I::Float64, E::Float64, numevals::Int64, x::Vector{Float64}, w::Vector{Float64}, gw::Vector{Float64}, n::Int64, atol::Float64, rtol::Float64, maxevals::Int64, nrm::typeof(LinearAlgebra.norm))
@ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:48
[3] do_quadgk(f::QuadGK.var"#14#23"{var"#1#2", Float64, Float64}, s::Tuple{Float64, Float64}, n::Int64, atol::Nothing, rtol::Nothing, maxevals::Int64, nrm::typeof(LinearAlgebra.norm), segbuf::Nothing)
@ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:37
[4] #46
@ ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:219 [inlined]
[5] handle_infinities(workfunc::QuadGK.var"#46#47"{Nothing, Nothing, Int64, Int64, typeof(LinearAlgebra.norm), Nothing}, f::var"#1#2", s::Tuple{Float64, Float64})
@ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:110
[6] quadgk(::Function, ::Float64, ::Vararg{Float64}; atol::Nothing, rtol::Nothing, maxevals::Int64, order::Int64, norm::Function, segbuf::Nothing)
@ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:218
[7] quadgk(::Function, ::Float64, ::Vararg{Float64})
@ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:216
[8] quadgk(::Function, ::Int64, ::Vararg{Any}; kws::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:213
[9] quadgk(::Function, ::Int64, ::Vararg{Any})
@ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:213
[10] top-level scope
@ In[2]:1
I think that the error is due to a call by quadgk of the integrand at the upper bound, which should not happen given the documentation of quadgk. https://docs.juliahub.com/QuadGK/hq5Ol/2.8.0/quadgk-examples/
This is confirmed by doing the change of variable done by quadgk manually (x=1+t/(1-t), dx = 1/(1-t)^2dt) as can be seen in the following code.
using QuadGK
function integrand(t)
if t==1.0
println(t)
end
return 1/(1+t/(1-t))^1.1 * 1/(1-t)^2
end
quadgk(integrand,0,1)
1.0
DomainError with 0.9999999999999964:
integrand produced NaN in the interval (0.9999999999999929, 1.0)
Stacktrace:
[1] evalrule(f::typeof(integrand), a::Float64, b::Float64, x::Vector{Float64}, w::Vector{Float64}, gw::Vector{Float64}, nrm::typeof(LinearAlgebra.norm))
@ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/evalrule.jl:37
[2] adapt(f::typeof(integrand), segs::Vector{QuadGK.Segment{Float64, Float64, Float64}}, I::Float64, E::Float64, numevals::Int64, x::Vector{Float64}, w::Vector{Float64}, gw::Vector{Float64}, n::Int64, atol::Float64, rtol::Float64, maxevals::Int64, nrm::typeof(LinearAlgebra.norm))
@ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:48
[3] do_quadgk(f::typeof(integrand), s::Tuple{Int64, Int64}, n::Int64, atol::Nothing, rtol::Nothing, maxevals::Int64, nrm::typeof(LinearAlgebra.norm), segbuf::Nothing)
@ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:37
[4] #46
@ ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:219 [inlined]
[5] handle_infinities(workfunc::QuadGK.var"#46#47"{Nothing, Nothing, Int64, Int64, typeof(LinearAlgebra.norm), Nothing}, f::typeof(integrand), s::Tuple{Int64, Int64})
@ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:118
[6] quadgk(::Function, ::Int64, ::Vararg{Int64}; atol::Nothing, rtol::Nothing, maxevals::Int64, order::Int64, norm::Function, segbuf::Nothing)
@ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:218
[7] quadgk(::Function, ::Int64, ::Vararg{Int64})
@ QuadGK ~/.julia/packages/QuadGK/XqIlh/src/adapt.jl:216
[8] top-level scope
@ In[15]:1
Reported by @sosock in PainterQubits/Unitful.jl#597 (comment):
julia> using Unitful, QuadGK
julia> quadgk(x->exp(-x/(1.0u"m")), 0.0u"m", Inf*u"m")
ERROR: DomainError with Inf m:
integrand produced NaN m in the interval (0.0 m, Inf m)
Stacktrace:
[...]
I think I might have discovered a recent regression in QuadGK
when using the quadgk
function with ranges specified in Unitful
types. I tried to dig into the source by following the stack trace, and it seems like it might be related to the recent changes to handle_infinities
, but I'm not familiar enough with this package or Julia's underlying type mechanics to nail down the issue at the moment.
I was able to reproduce the issue using the following minimal code snippet. This snippet should integrate a constant 1 Joule value over the domain [0,10] seconds, which should result in a value of 10 Joule-seconds.
using QuadGK, Unitful
f(t) = 1.0u"J"
quadgk(f, 0.0u"s", 10.0u"s")
System information for this configuration:
versioninfo()
# Julia Version 1.8.5
# Commit 17cfb8e65e (2023-01-08 06:45 UTC)
# Platform Info:
# OS: Windows (x86_64-w64-mingw32)
# CPU: 32 × AMD Ryzen 9 3950X 16-Core Processor
# WORD_SIZE: 64
# LIBM: libopenlibm
# LLVM: libLLVM-13.0.1 (ORCJIT, znver2)
# Threads: 1 on 32 virtual cores
import Pkg
Pkg.status(["Unitful", "QuadGK"])
# Status `C:\Users\mikei\.julia\environments\v1.8\Project.toml`
# [1fd47b50] QuadGK v2.6.0
# [1986cc42] Unitful v1.12.2
Running the code snippet above produces an error
using QuadGK, Unitful
f(t) = 1.0u"J"
# f (generic function with 1 method)
quadgk(f, 0.0u"s", 10.0u"s")
#=
ERROR: MethodError: no method matching Quantity{Float64, 𝐓 , Unitful.FreeUnits{(s,), 𝐓 , nothing}}(::Int64)
The applicable method may be too new: running in world age 32472, while current world is 32475.
Closest candidates are:
Quantity{T, D, U}(::Number) where {T, D, U} at C:\Users\mikei\.julia\packages\Unitful\fbiNW\src\types.jl:151 (method too new to be called from this world context.)
(::Type{T})(::T) where T<:Number at boot.jl:772
(::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
...
Stacktrace:
[1] convert(#unused#::Type{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}} , x::Int64)
@ Base .\number.jl:7
[2] one(#unused#::Type{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}} )
@ Base .\number.jl:334
[3] #s14#1
@ C:\Users\mikei\.julia\packages\QuadGK\kf0xA\src\gausskronrod.jl:260 [inlined]
[4] var"#s14#1"(T::Any, ::Any, #unused#::Any, n::Any)
@ QuadGK .\none:0
[5] (::Core.GeneratedFunctionStub)(::Any, ::Vararg{Any})
@ Core .\boot.jl:582
[6] do_quadgk(f::typeof(f), s::Tuple{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}, Quantity{Float64, 𝐓
, Unitful.FreeUnits{(s,), 𝐓, nothing}}} , n::Int64, atol::Nothing, rtol::Nothing, maxevals::Int64, nrm::typeof(LinearAlgebra.norm), segbuf::Nothing)
@ QuadGK C:\Users\mikei\.julia\packages\QuadGK\kf0xA\src\adapt.jl:7
[7] #28
@ C:\Users\mikei\.julia\packages\QuadGK\kf0xA\src\adapt.jl:186 [inlined]
[8] handle_infinities(workfunc::QuadGK.var"#28#29"{Nothing, Nothing, Int64, Int64, typeof(LinearAlgebra.norm), Nothing}, f::typeof(f), s::Tuple{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}, Quantity{Float64, 𝐓, Unitful.Fre
eUnits{(s,), 𝐓, nothing}}} )
@ QuadGK C:\Users\mikei\.julia\packages\QuadGK\kf0xA\src\adapt.jl:115
[9] quadgk(::Function, ::Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}} , ::Vararg{Quantity{Float64, 𝐓, U
nitful.FreeUnits{(s,), 𝐓, nothing}}} ; atol::Nothing, rtol::Nothing, maxevals::Int64, order::Int64, norm::Function, segbuf::Nothing)
@ QuadGK C:\Users\mikei\.julia\packages\QuadGK\kf0xA\src\adapt.jl:185
[10] quadgk(::Function, ::Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}} , ::Vararg{Quantity{Float64, 𝐓, U
nitful.FreeUnits{(s,), 𝐓, nothing}}} )
@ QuadGK C:\Users\mikei\.julia\packages\QuadGK\kf0xA\src\adapt.jl:183
[11] top-level scope
@ REPL[3]:1
=#
For the next configuration, I simply pinned QuadGK
to version 2.5.0
in Julia's Pkg
REPL and re-ran the test. System information for this configuration:
versioninfo()
# Julia Version 1.8.5
# Commit 17cfb8e65e (2023-01-08 06:45 UTC)
# Platform Info:
# OS: Windows (x86_64-w64-mingw32)
# CPU: 32 × AMD Ryzen 9 3950X 16-Core Processor
# WORD_SIZE: 64
# LIBM: libopenlibm
# LLVM: libLLVM-13.0.1 (ORCJIT, znver2)
# Threads: 1 on 32 virtual cores
import Pkg
Pkg.status(["Unitful", "QuadGK"])
# ⌃ [1fd47b50] QuadGK v2.5.0 ⚲
# [1986cc42] Unitful v1.12.2
# Info Packages marked with ⌃ have new versions available and may be upgradable.
In this configuration, the test runs as expected.
using QuadGK, Unitful
f(t) = 1.0u"J"
# f (generic function with 1 method)
quadgk(f, 0.0u"s", 10.0u"s")
# (10.0 J s, 0.0 J s)
I thought this was supported, but then this integral should give 2πi by the residue theorem:
julia> quadgk(z -> cos(z)/z, 1, im, -1, -im)
(0.50046301807924 + 4.71238898038469im, 4.027482497221434e-9)
and the result looks clearly wrong. 😢
In the documentation for kronrod()
I see references to both wg
and gw
. Also there are two entries for kronrod()
on that page, the second near the bottom.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.