1 module MultiQuad  2 3 using QuadGK, Cuba, Unitful  4 5 export quad, dblquad, tplquad  6 7 @doc raw"""  8  quad(arg::Function, x1, x2; method = :quadgk, kwargs...)  9 10 Performs the integral \int_{x1}^{x2}f(x)dx  11 See [QuadGK](https://github.com/JuliaMath/QuadGK.jl) and [Cuba.jl](https://giordano.github.io/Cuba.jl/stable/) for all the available keyword arguments.  12 13 # Examples  14 jldoctest  15 func(x) = x^2*exp(-x)  16 integral, error = quad(func, 0, 4)  17   18 19 quad can handle units of measurement through the [Unitful](https://github.com/PainterQubits/Unitful.jl) package:  20 21 jldoctest  22 using Unitful  23 24 func(x) = x^2  25 integral, error = quad(func, 1u"m", 5u"m")  26   27 """  28 function quad(arg::Function, x1, x2; method = :quadgk, kwargs...)  29 30 24  if method == :suave  31 24  integrate = suave  32 24  elseif method == :vegas  33 24  integrate = vegas  34 24  elseif method == :quadgk  35 24  integrate = quadgk  36  else  37 18  ex = ErrorException("Integration method $method is not supported!")  38 24  throw(ex)  39  end  40 41 24  if method == :quadgk  42 24  return quadgk(arg, x1, x2; kwargs...)  43  end  44 45 18  units = unit(arg(x1)) * unit(x1)  46 47 24  arg2(a) = ustrip(units, (x2 - x1) * arg((x2 - x1) * a + x1))::Float64  48 49 24  function integrand(x, f)  50 24  f[1] = arg2(x[1])  51  end  52 53 24  result, err = integrate(integrand, 1, 1; kwargs...)  54 55 6  if units == Unitful.NoUnits  56 24  return result[1], err[1]  57  else  58 24  return (result[1], err[1]) .* units  59  end  60 end  61 62 @doc raw"""  63  dblquad(arg::Function, x1, x2, y1::Function, y2::Function; method = :cuhre, kwargs...)  64 65 Performs the integral \int_{x1}^{x2}\int_{y1(x)}^{y2(x)}f(y,x)dydx  66 See [Cuba.jl](https://giordano.github.io/Cuba.jl/stable/) for all the available keyword arguments.  67 68 # Examples  69 jldoctest  70 func(y,x) = sin(x)*y^2  71 integral, error = dblquad(func, 1, 2, x->0, x->x^2, rtol=1e-9)  72   73 74 dblquad can handle units of measurement through the [Unitful](https://github.com/PainterQubits/Unitful.jl) package:  75 76 jldoctest  77 using Unitful  78 79 func(y,x) = x*y^2  80 integral, error = dblquad(func, 1u"m", 2u"m", x->0u"m^2", x->x^2, rtol=1e-9)  81   82 """  83 function dblquad(  84  arg::Function,  85  x1,  86  x2,  87  y1::Function,  88  y2::Function;  89  method = :cuhre,  90  kwargs...,  91 )  92 93 24  if method == :cuhre  94 24  integrate = cuhre  95 24  elseif method == :divonne  96 24  integrate = divonne  97 24  elseif method == :suave  98 24  integrate = suave  99 24  elseif method == :vegas  100 24  integrate = vegas  101  else  102 18  ex = ErrorException("Integration method$method is not supported!")  103 24  throw(ex)  104  end  105 106 24  units = unit(arg(y1(x1), x1)) * unit(x1) * unit(y1(x1))  107 108 24  arg1(a, x) = (y2(x) - y1(x)) * arg((y2(x) - y1(x)) * a + y1(x), x)  109 110 111 24  arg2(a, b) = ustrip(units, (x2 - x1) * arg1(a, (x2 - x1) * b + x1))::Float64  112 113 24  function integrand(x, f)  114 24  f[1] = arg2(x[1], x[2])  115  end  116 117 24  result, err = integrate(integrand, 2, 1; kwargs...)  118 119 6  if units == Unitful.NoUnits  120 24  return result[1], err[1]  121  else  122 24  return (result[1], err[1]) .* units  123  end  124 end  125 126 @doc raw"""  127  tplquad(arg::Function, x1, x2, y1::Function, y2::Function, z1::Function, z2::Function; method = :cuhre, kwargs...)  128 129 Performs the integral \int_{x1}^{x2}\int_{y1(x)}^{y2(x)}\int_{z1(x,y)}^{z2(x,y)}f(z,y,x)dzdydx  130 See [Cuba.jl](https://giordano.github.io/Cuba.jl/stable/) for all the available keyword arguments.  131 132 # Examples  133 jldoctest  134 func(z,y,x) = sin(z)*y*x  135 integral, error = tplquad(func, 0, 4, x->x, x->x^2, (x,y)->2, (x,y)->3*x)  136   137 138 tplquad can handle units of measurement through the [Unitful](https://github.com/PainterQubits/Unitful.jl) package:  139 140 jldoctest  141 using Unitful  142 143 func(z,y,x) = sin(z)*y*x  144 integral, error = tplquad(func, 0u"m", 4u"m", x->0u"m^2", x->x^2, (x,y)->0, (x,y)->3)  145   146 """  147 function tplquad(  148  arg::Function,  149  x1,  150  x2,  151  y1::Function,  152  y2::Function,  153  z1::Function,  154  z2::Function;  155  method = :cuhre,  156  kwargs...,  157 )  158 159 24  if method == :cuhre  160 24  integrate = cuhre  161 24  elseif method == :divonne  162 24  integrate = divonne  163 24  elseif method == :suave  164 24  integrate = suave  165 24  elseif method == :vegas  166 24  integrate = vegas  167  else  168 0  ex = ErrorException("Integration method \$method is not supported!")  169 0  throw(ex)  170  end  171 172 24  units = unit(arg(z1(x1, y1(x1)), y1(x1), x1)) * unit(x1) * unit(y1(x1)) *  173  unit(z1(y1(x1), x1))  174 175 24  arg0(a, y, x) =  176  (z2(x, y) - z1(x, y)) * arg((z2(x, y) - z1(x, y)) * a + z1(x, y), y, x)  177 178 24  arg1(a, b, x) = (y2(x) - y1(x)) * arg0(a, (y2(x) - y1(x)) * b + y1(x), x)  179 180 181 24  arg2(a, b, c) =  182  ustrip(units, (x2 - x1) * arg1(a, b, (x2 - x1) * c + x1))::Float64  183 184 24  function integrand(x, f)  185 24  f[1] = arg2(x[1], x[2], x[3])  186  end  187 188 24  result, err = integrate(integrand, 3, 1; kwargs...)  189 190 6  if units == Unitful.NoUnits  191 24  return result[1], err[1]  192  else  193 24  return (result[1], err[1]) .* units  194  end  195 end  196 197 198 end # module 

Read our documentation on viewing source code .