aurelio-amerio / MultiQuad.jl
Showing 1 of 7 files from the diff.
Other files ignored by Codecov
LICENSE has changed.
test.jl is new.
test/runtests.jl has changed.
Project.toml has changed.
.travis.yml was deleted.
README.md has changed.

@@ -1,631 +1,677 @@
Loading
1 -
module MultiQuad
2 -
3 -
using QuadGK, Cuba, HCubature, Unitful, FastGaussQuadrature, LinearAlgebra, LRUCache
4 -
5 -
export quad, dblquad, tplquad
6 -
7 -
8 -
function _quadgk_wrapper_1D(arg::Function, x1, x2; kwargs...)
9 -
    return quadgk(arg, x1, x2; kwargs...)
10 -
end
11 -
12 -
function _cuba_wrapper_1D(arg::Function, x1, x2; method::Symbol, kwargs...)
13 -
    if method == :suave
14 -
        integrate = suave
15 -
    elseif method == :vegas
16 -
        integrate = vegas
17 -
    else
18 -
        ex = ErrorException("Integration method $method is not supported!")
19 -
        throw(ex)
20 -
    end
21 -
22 -
    units = unit(arg(x1)) * unit(x1)
23 -
24 -
    arg2(a) = ustrip(units, (x2 - x1) * arg((x2 - x1) * a + x1))::Float64
25 -
26 -
    function integrand(x, f)
27 -
        f[1] = arg2(x[1])
28 -
    end
29 -
30 -
    result, err = integrate(integrand, 1, 1; kwargs...)
31 -
32 -
    if units == Unitful.NoUnits
33 -
        return result[1], err[1]
34 -
    else
35 -
        return (result[1], err[1]) .* units
36 -
    end
37 -
end
38 -
39 -
#order, method, optional
40 -
41 -
function get_weights_fastgaussquadrature(order::Int, method::Symbol, optional::Real=NaN)
42 -
    if method == :gausslegendre
43 -
        xw = gausslegendre
44 -
        return xw(order)
45 -
46 -
    elseif method == :gausshermite
47 -
        xw = gausshermite
48 -
        return xw(order)
49 -
50 -
    elseif method == :gausslaguerre
51 -
        xw = gausslaguerre
52 -
        if isnan(optional)
53 -
            return xw(order)
54 -
        else
55 -
            α = optional
56 -
            return xw(order, α)
57 -
        end
58 -
59 -
    elseif method == :gausschebyshev
60 -
        xw = gausschebyshev
61 -
        if isnan(optional)
62 -
            return xw(order)
63 -
        else
64 -
            kind = optional
65 -
            return xw(order, kind)
66 -
        end
67 -
68 -
    # elseif method == :gaussjacobi
69 -
    #     xw = gaussjacobi
70 -
    #     return xw(order)
71 -
72 -
    elseif method == :gaussradau
73 -
        xw = gaussradau
74 -
        return xw(order)
75 -
76 -
    elseif method == :gausslobatto
77 -
        xw = gausslobatto
78 -
        return xw(order)
79 -
    
80 -
    else
81 -
        ex = ErrorException("Integration method $method is not supported!")
82 -
        throw(ex)
83 -
    end
84 -
85 -
    
86 -
end
87 -
88 -
const lru_xw_fastgaussquadrature = LRU{
89 -
    Tuple{Int,Symbol,Real},
90 -
    Tuple{Vector{Float64},Vector{Float64}},
91 -
}(maxsize = Int(1e5))
92 -
93 -
function get_weights_fastgaussquadrature_cached(order::Int, method::Symbol, optional::Real=NaN)
94 -
    get!(lru_xw_fastgaussquadrature, (order, method, optional)) do
95 -
        get_weights_fastgaussquadrature(order, method, optional)
96 -
    end
97 -
end
98 -
99 -
function _fastgaussquadrature_wrapper_1D(arg::Function, x1, x2; method::Symbol, order::Int, kwargs...)
100 -
    if ! (method in [:gausslegendre, :gausschebyshev, :gaussradau, :gausslobatto])
101 -
        ex = ErrorException("Integration method $method is not supported!")
102 -
        throw(ex)
103 -
    end
104 -
105 -
    if haskey(kwargs, :α)
106 -
        optional = kwargs[:α]
107 -
    elseif haskey(kwargs, :kind)
108 -
        optional = kwargs[:kind]
109 -
    else 
110 -
        optional = NaN
111 -
    end
112 -
113 -
    x, w = get_weights_fastgaussquadrature_cached(order, method, optional)
114 -
115 -
    b=x2
116 -
    a=x1
117 -
118 -
    b_a_2 = (b-a)/2
119 -
    a_plus_b_2 = (a+b)/2
120 -
121 -
    arg_gauss(x) = arg(b_a_2*x + a_plus_b_2)
122 -
    
123 -
    integral = b_a_2 * dot(w, arg_gauss.(x))
124 -
    return integral, NaN
125 -
126 -
end
127 -
128 -
function _fastgaussquadrature_wrapper_1D(arg::Function; method::Symbol, order::Int, kwargs...)
129 -
    if ! (method in [:gausshermite , :gausslaguerre])
130 -
        ex = ErrorException("Integration method $method is not supported!")
131 -
        throw(ex)
132 -
    end
133 -
134 -
    if haskey(kwargs, :α)
135 -
        optional = kwargs[:α]
136 -
    elseif haskey(kwargs, :kind)
137 -
        optional = kwargs[:kind]
138 -
    else 
139 -
        optional = NaN
140 -
    end
141 -
142 -
    x, w = get_weights_fastgaussquadrature_cached(order, method, optional)
143 -
    
144 -
    integral = dot(w, arg.(x))
145 -
    return integral, NaN
146 -
147 -
end
148 -
149 -
@doc raw"""
150 -
    quad(arg::Function, x1, x2[; method = :quadgk, oder=1000, kwargs...])
151 -
152 -
Performs the integral ``\int_{x1}^{x2}f(x)dx``
153 -
154 -
Available integration methods:
155 -
- `:suave`
156 -
- `:vegas`
157 -
- `:quadgk`
158 -
159 -
There are several fixed order quadrature methods available based on [`FastGaussQuadrature`](https://github.com/JuliaApproximation/FastGaussQuadrature.jl)
160 -
161 -
- `:gausslegendre` 
162 -
- `:gausshermite` 
163 -
- `:gausslaguerre` 
164 -
- `:gausschebyshev` 
165 -
- `:gaussradau` 
166 -
- `:gausslobatto`
167 -
168 -
If a specific integration routine is not needed, it is suggested to use `:quadgk` (the default option).
169 -
For highly oscillating integrals, it is advised to use `:gausslegendre`with a high order (~10000).
170 -
171 -
Please note that `:gausshermite` and `:gausslaguerre` are used to specific kind of integrals with infinite bounds.
172 -
173 -
For more informations on these integration techniques, please follow the official documentation [here](https://juliaapproximation.github.io/FastGaussQuadrature.jl/stable/gaussquadrature/).
174 -
175 -
See [QuadGK](https://github.com/JuliaMath/QuadGK.jl) and [Cuba.jl](https://giordano.github.io/Cuba.jl/stable/) for all the available keyword arguments.
176 -
# Examples
177 -
```jldoctest
178 -
using MultiQuad
179 -
180 -
func(x) = x^2*exp(-x)
181 -
quad(func, 0, 4) # equivalent to quad(func, 0, 4, method=:quadgk)
182 -
quad(func, 0, 4, method=:gausslegendre, order=1000)
183 -
184 -
# for certain kinds of integrals with infinite bounds, it may be possible to use a specific integration routine
185 -
func(x) = x^2*exp(-x)
186 -
quad(func, 0, Inf, method=:quadgk)[1]
187 -
# gausslaguerre computes the integral of f(x)*exp(-x) from 0 to infinity
188 -
quad(x -> x^4, method=:gausslaguerre, order=10000)[1] 
189 -
```
190 -
191 -
`quad` can handle units of measurement through the [`Unitful`](https://github.com/PainterQubits/Unitful.jl) package:
192 -
193 -
```jldoctest
194 -
using Unitful
195 -
196 -
func(x) = x^2
197 -
integral, error = quad(func, 1u"m", 5u"m")
198 -
```
199 -
"""
200 -
function quad(arg::Function, x1, x2; method::Symbol = :quadgk, order::Int=1000, kwargs...)
201 -
202 -
    if method in [:suave, :vegas]
203 -
        return _cuba_wrapper_1D(arg, x1, x2; method=method)
204 -
    elseif method == :quadgk
205 -
        return _quadgk_wrapper_1D(arg, x1, x2; kwargs...)
206 -
    elseif method in [:gausslegendre, :gausschebyshev, :gaussradau, :gausslobatto]
207 -
        return _fastgaussquadrature_wrapper_1D(arg, x1, x2; method=method, order=order, kwargs...)    
208 -
    else
209 -
        ex = ErrorException("Integration method $method is not supported!")
210 -
        throw(ex)
211 -
    end
212 -
end
213 -
214 -
@doc raw"""
215 -
    quad(arg::Function; method[, oder=1000, kwargs...])
216 -
"""
217 -
function quad(arg::Function; method::Symbol, order::Int=1000, kwargs...)
218 -
219 -
    if method in [:gausshermite , :gausslaguerre]
220 -
        return _fastgaussquadrature_wrapper_1D(arg; method=method, order=order, kwargs...)    
221 -
    else
222 -
        ex = ErrorException("Integration method $method is not supported!")
223 -
        throw(ex)
224 -
    end
225 -
end
226 -
227 -
function _cuba_wrapper_2D(arg::Function, x1, x2, y1::Function, y2::Function; method::Symbol, kwargs...)
228 -
    if method == :cuhre
229 -
        integrate = cuhre
230 -
    elseif method == :divonne
231 -
        integrate = divonne
232 -
    elseif method == :suave
233 -
        integrate = suave
234 -
    elseif method == :vegas
235 -
        integrate = vegas
236 -
    else
237 -
        ex = ErrorException("Integration method $method is not supported!")
238 -
        throw(ex)
239 -
    end
240 -
241 -
    units = unit(arg(y1(x1), x1)) * unit(x1) * unit(y1(x1))
242 -
243 -
    arg1(a, x) = (y2(x) - y1(x)) * arg((y2(x) - y1(x)) * a + y1(x), x)
244 -
245 -
246 -
    arg2(a, b) = ustrip(units, (x2 - x1) * arg1(a, (x2 - x1) * b + x1))::Float64
247 -
248 -
    function integrand2(x, f)
249 -
        f[1] = arg2(x[1], x[2])
250 -
    end
251 -
252 -
    result, err = integrate(integrand2, 2, 1; kwargs...)
253 -
254 -
    if units == Unitful.NoUnits
255 -
        return result[1], err[1]
256 -
    else
257 -
        return (result[1], err[1]) .* units
258 -
    end
259 -
end
260 -
261 -
function _cuba_wrapper_2D(arg::Function, x1, x2, y1, y2; method::Symbol, kwargs...)
262 -
    if method == :cuhre
263 -
        integrate = cuhre
264 -
    elseif method == :divonne
265 -
        integrate = divonne
266 -
    elseif method == :suave
267 -
        integrate = suave
268 -
    elseif method == :vegas
269 -
        integrate = vegas
270 -
    else
271 -
        ex = ErrorException("Integration method $method is not supported!")
272 -
        throw(ex)
273 -
    end
274 -
275 -
    units = unit(arg(y1, x1)) * unit(x1) * unit(y1)
276 -
277 -
    arg1(a, x) = (y2 - y1) * arg((y2 - y1) * a + y1, x)
278 -
279 -
    arg2(a, b) = ustrip(units, (x2 - x1) * arg1(a, (x2 - x1) * b + x1))::Float64
280 -
281 -
    function integrand2(x, f)
282 -
        f[1] = arg2(x[1], x[2])
283 -
    end
284 -
285 -
    result, err = integrate(integrand2, 2, 1; kwargs...)
286 -
287 -
    if units == Unitful.NoUnits
288 -
        return result[1], err[1]
289 -
    else
290 -
        return (result[1], err[1]) .* units
291 -
    end
292 -
end
293 -
294 -
295 -
function _hcubature_wrapper_2D(arg::Function, x1, x2, y1::Function, y2::Function; kwargs...)
296 -
    units = unit(arg(y1(x1), x1)) * unit(x1) * unit(y1(x1))
297 -
298 -
    arg1(a, x) = (y2(x) - y1(x)) * arg((y2(x) - y1(x)) * a + y1(x), x)
299 -
300 -
301 -
    arg2(a, b) = ustrip(units, (x2 - x1) * arg1(a, (x2 - x1) * b + x1))::Float64
302 -
303 -
304 -
    function integrand(arr)
305 -
        return arg2(arr[1], arr[2])
306 -
    end
307 -
308 -
    min_arr = [0, 0]
309 -
    max_arr = [1, 1]
310 -
    result, err = hcubature(integrand, min_arr, max_arr; kwargs...)
311 -
    
312 -
313 -
    if units == Unitful.NoUnits
314 -
        return result[1], err[1]
315 -
    else
316 -
        return (result[1], err[1]) .* units
317 -
    end
318 -
end
319 -
320 -
function _hcubature_wrapper_2D(arg::Function, x1, x2, y1, y2; kwargs...)
321 -
    units = unit(arg(y1, x1)) * unit(x1) * unit(y1)
322 -
323 -
    arg1(a, x) = (y2 - y1) * arg((y2 - y1) * a + y1, x)
324 -
325 -
    arg2(a, b) = ustrip(units, (x2 - x1) * arg1(a, (x2 - x1) * b + x1))::Float64
326 -
327 -
328 -
    function integrand(arr)
329 -
        return arg2(arr[1], arr[2])
330 -
    end
331 -
332 -
    min_arr = [0, 0]
333 -
    max_arr = [1, 1]
334 -
    result, err = hcubature(integrand, min_arr, max_arr; kwargs...)
335 -
    
336 -
337 -
    if units == Unitful.NoUnits
338 -
        return result[1], err[1]
339 -
    else
340 -
        return (result[1], err[1]) .* units
341 -
    end
342 -
end
343 -
344 -
@doc raw"""
345 -
    dblquad(arg::Function, x1, x2, y1::Function, y2::Function; method = :cuhre, kwargs...)
346 -
347 -
Performs the integral ``\int_{x1}^{x2}\int_{y1(x)}^{y2(x)}f(y,x)dydx``
348 -
349 -
Available integration methods:
350 -
- `:cuhre`
351 -
- `:divonne`
352 -
- `:suave`
353 -
- `:vegas`
354 -
- `:hcubature`
355 -
356 -
See [Cuba.jl](https://giordano.github.io/Cuba.jl/stable/) for all the available keyword arguments fro the `:cuhre`, `:divonne`, `:suave` and `:vegas` methods.
357 -
See [HCubature](https://github.com/stevengj/HCubature.jl) for all the available keywords for the `:hcubature` method.
358 -
359 -
# Examples
360 -
```jldoctest
361 -
func(y,x) = sin(x)*y^2
362 -
integral, error = dblquad(func, 1, 2, x->0, x->x^2, rtol=1e-9)
363 -
```
364 -
365 -
`dblquad` can handle units of measurement through the [`Unitful`](https://github.com/PainterQubits/Unitful.jl) package:
366 -
367 -
```jldoctest
368 -
using Unitful
369 -
370 -
func(y,x) = x*y^2
371 -
integral, error = dblquad(func, 1u"m", 2u"m", x->0u"m^2", x->x^2, rtol=1e-9)
372 -
```
373 -
"""
374 -
function dblquad(
375 -
    arg::Function,
376 -
    x1,
377 -
    x2,
378 -
    y1::Function,
379 -
    y2::Function;
380 -
    method = :hcubature,
381 -
    kwargs...,
382 -
)
383 -
384 -
    if method in [:cuhre, :divonne, :suave, :vegas]
385 -
        return _cuba_wrapper_2D(arg, x1, x2, y1, y2; method=method)
386 -
    elseif method == :hcubature
387 -
        return _hcubature_wrapper_2D(arg, x1, x2, y1, y2; kwargs...)
388 -
    else
389 -
        ex = ErrorException("Integration method $method is not supported!")
390 -
        throw(ex)
391 -
    end
392 -
end
393 -
394 -
function dblquad(
395 -
    arg::Function,
396 -
    x1,
397 -
    x2,
398 -
    y1,
399 -
    y2;
400 -
    method = :hcubature,
401 -
    kwargs...,
402 -
)
403 -
404 -
    if method in [:cuhre, :divonne, :suave, :vegas]
405 -
        return _cuba_wrapper_2D(arg, x1, x2, y1, y2; method=method)
406 -
    elseif method == :hcubature
407 -
        return _hcubature_wrapper_2D(arg, x1, x2, y1, y2; kwargs...)
408 -
    else
409 -
        ex = ErrorException("Integration method $method is not supported!")
410 -
        throw(ex)
411 -
    end
412 -
end
413 -
414 -
function _cuba_wrapper_3D(arg::Function, x1, x2, y1::Function, y2::Function, z1::Function, z2::Function;  method::Symbol, kwargs...)
415 -
    if method == :cuhre
416 -
        integrate = cuhre
417 -
    elseif method == :divonne
418 -
        integrate = divonne
419 -
    elseif method == :suave
420 -
        integrate = suave
421 -
    elseif method == :vegas
422 -
        integrate = vegas
423 -
    else
424 -
        ex = ErrorException("Integration method $method is not supported!")
425 -
        throw(ex)
426 -
    end
427 -
428 -
    units = unit(arg(z1(x1, y1(x1)), y1(x1), x1)) * unit(x1) * unit(y1(x1)) *
429 -
            unit(z1(y1(x1), x1))
430 -
431 -
    arg0(a, y, x) =
432 -
        (z2(x, y) - z1(x, y)) * arg((z2(x, y) - z1(x, y)) * a + z1(x, y), y, x)
433 -
434 -
    arg1(a, b, x) = (y2(x) - y1(x)) * arg0(a, (y2(x) - y1(x)) * b + y1(x), x)
435 -
436 -
437 -
    arg2(a, b, c) =
438 -
        ustrip(units, (x2 - x1) * arg1(a, b, (x2 - x1) * c + x1))::Float64
439 -
440 -
    
441 -
    function integrand2(x, f)
442 -
        f[1] = arg2(x[1], x[2], x[3])
443 -
    end
444 -
445 -
    result, err = integrate(integrand2, 3, 1; kwargs...)
446 -
    
447 -
448 -
    if units == Unitful.NoUnits
449 -
        return result[1], err[1]
450 -
    else
451 -
        return (result[1], err[1]) .* units
452 -
    end
453 -
end
454 -
455 -
function _cuba_wrapper_3D(arg::Function, x1, x2, y1, y2, z1, z2;  method::Symbol, kwargs...)
456 -
    if method == :cuhre
457 -
        integrate = cuhre
458 -
    elseif method == :divonne
459 -
        integrate = divonne
460 -
    elseif method == :suave
461 -
        integrate = suave
462 -
    elseif method == :vegas
463 -
        integrate = vegas
464 -
    else
465 -
        ex = ErrorException("Integration method $method is not supported!")
466 -
        throw(ex)
467 -
    end
468 -
469 -
    units = unit(arg(z1, y1, x1)) * unit(x1) * unit(y1) * unit(z1)
470 -
471 -
    arg0(a, y, x) = (z2 - z1) * arg((z2 - z1) * a + z1, y, x)
472 -
473 -
    arg1(a, b, x) = (y2 - y1) * arg0(a, (y2 - y1) * b + y1, x)
474 -
475 -
476 -
    arg2(a, b, c) =
477 -
        ustrip(units, (x2 - x1) * arg1(a, b, (x2 - x1) * c + x1))::Float64
478 -
479 -
    
480 -
    function integrand2(x, f)
481 -
        f[1] = arg2(x[1], x[2], x[3])
482 -
    end
483 -
484 -
    result, err = integrate(integrand2, 3, 1; kwargs...)
485 -
    
486 -
487 -
    if units == Unitful.NoUnits
488 -
        return result[1], err[1]
489 -
    else
490 -
        return (result[1], err[1]) .* units
491 -
    end
492 -
end
493 -
494 -
function _hcubature_wrapper_3D(arg::Function, x1, x2, y1::Function, y2::Function, z1::Function, z2::Function; kwargs...)
495 -
496 -
    units = unit(arg(z1(x1, y1(x1)), y1(x1), x1)) * unit(x1) * unit(y1(x1)) *
497 -
            unit(z1(y1(x1), x1))
498 -
499 -
    arg0(a, y, x) =
500 -
        (z2(x, y) - z1(x, y)) * arg((z2(x, y) - z1(x, y)) * a + z1(x, y), y, x)
501 -
502 -
    arg1(a, b, x) = (y2(x) - y1(x)) * arg0(a, (y2(x) - y1(x)) * b + y1(x), x)
503 -
504 -
505 -
    arg2(a, b, c) =
506 -
        ustrip(units, (x2 - x1) * arg1(a, b, (x2 - x1) * c + x1))::Float64
507 -
508 -
    
509 -
    function integrand(arr)
510 -
        return arg2(arr[1], arr[2], arr[3])
511 -
    end
512 -
513 -
    min_arr = [0, 0, 0]
514 -
    max_arr = [1, 1, 1]
515 -
    result, err = hcubature(integrand, min_arr, max_arr; kwargs...)
516 -
    
517 -
    if units == Unitful.NoUnits
518 -
        return result[1], err[1]
519 -
    else
520 -
        return (result[1], err[1]) .* units
521 -
    end
522 -
end
523 -
524 -
function _hcubature_wrapper_3D(arg::Function, x1, x2, y1, y2, z1, z2; kwargs...)
525 -
526 -
    units = unit(arg(z1, y1, x1)) * unit(x1) * unit(y1) * unit(z1)
527 -
528 -
    arg0(a, y, x) = (z2 - z1) * arg((z2 - z1) * a + z1, y, x)
529 -
530 -
    arg1(a, b, x) = (y2 - y1) * arg0(a, (y2 - y1) * b + y1, x)
531 -
532 -
533 -
    arg2(a, b, c) =
534 -
        ustrip(units, (x2 - x1) * arg1(a, b, (x2 - x1) * c + x1))::Float64
535 -
536 -
    
537 -
    function integrand(arr)
538 -
        return arg2(arr[1], arr[2], arr[3])
539 -
    end
540 -
541 -
    min_arr = [0, 0, 0]
542 -
    max_arr = [1, 1, 1]
543 -
    result, err = hcubature(integrand, min_arr, max_arr; kwargs...)
544 -
    
545 -
    if units == Unitful.NoUnits
546 -
        return result[1], err[1]
547 -
    else
548 -
        return (result[1], err[1]) .* units
549 -
    end
550 -
end
551 -
552 -
@doc raw"""
553 -
    tplquad(arg::Function, x1, x2, y1::Function, y2::Function, z1::Function, z2::Function; method = :cuhre, kwargs...)
554 -
555 -
Performs the integral ``\int_{x1}^{x2}\int_{y1(x)}^{y2(x)}\int_{z1(x,y)}^{z2(x,y)}f(z,y,x)dzdydx``
556 -
557 -
Available integration methods:
558 -
- `:cuhre`
559 -
- `:divonne`
560 -
- `:suave`
561 -
- `:vegas`
562 -
- `:hcubature`
563 -
564 -
See [Cuba.jl](https://giordano.github.io/Cuba.jl/stable/) for all the available keyword arguments fro the `:cuhre`, `:divonne`, `:suave` and `:vegas` methods.
565 -
See [HCubature](https://github.com/stevengj/HCubature.jl) for all the available keywords for the `:hcubature` method.
566 -
567 -
# Examples
568 -
```jldoctest
569 -
func(z,y,x) = sin(z)*y*x
570 -
integral, error = tplquad(func, 0, 4, x->x, x->x^2, (x,y)->2, (x,y)->3*x)
571 -
```
572 -
573 -
`tplquad` can handle units of measurement through the [`Unitful`](https://github.com/PainterQubits/Unitful.jl) package:
574 -
575 -
```jldoctest
576 -
using Unitful
577 -
578 -
func(z,y,x) = sin(z)*y*x
579 -
integral, error = tplquad(func, 0u"m", 4u"m", x->0u"m^2", x->x^2, (x,y)->0, (x,y)->3)
580 -
```
581 -
"""
582 -
function tplquad(
583 -
    arg::Function,
584 -
    x1,
585 -
    x2,
586 -
    y1::Function,
587 -
    y2::Function,
588 -
    z1::Function,
589 -
    z2::Function;
590 -
    method = :hcubature,
591 -
    kwargs...,
592 -
)
593 -
594 -
    if method in [:cuhre, :divonne, :suave, :vegas]
595 -
        return _cuba_wrapper_3D(arg, x1, x2, y1, y2, z1, z2;  method=method, kwargs...)
596 -
597 -
    elseif method == :hcubature
598 -
        return _hcubature_wrapper_3D(arg, x1, x2, y1, y2, z1, z2; kwargs...)
599 -
600 -
    else
601 -
        ex = ErrorException("Integration method $method is not supported!")
602 -
        throw(ex)
603 -
    end
604 -
end
605 -
606 -
function tplquad(
607 -
    arg::Function,
608 -
    x1,
609 -
    x2,
610 -
    y1,
611 -
    y2,
612 -
    z1,
613 -
    z2;
614 -
    method = :hcubature,
615 -
    kwargs...,
616 -
)
617 -
618 -
    if method in [:cuhre, :divonne, :suave, :vegas]
619 -
        return _cuba_wrapper_3D(arg, x1, x2, y1, y2, z1, z2;  method=method, kwargs...)
620 -
621 -
    elseif method == :hcubature
622 -
        return _hcubature_wrapper_3D(arg, x1, x2, y1, y2, z1, z2; kwargs...)
623 -
624 -
    else
625 -
        ex = ErrorException("Integration method $method is not supported!")
626 -
        throw(ex)
627 -
    end
628 -
end
629 -
630 -
631 -
end # module
1 +
module MultiQuad
2 +
3 +
using QuadGK, Cuba, HCubature, Unitful, FastGaussQuadrature, LinearAlgebra, LRUCache, Base.Threads, ProgressMeter
4 +
5 +
export quad, dblquad, tplquad
6 +
7 +
8 +
function _quadgk_wrapper_1D(arg::Function, x1, x2; kwargs...)
9 +
    return quadgk(arg, x1, x2; kwargs...)
10 +
end
11 +
12 +
function _cuba_wrapper_1D(arg::Function, x1, x2; method::Symbol, kwargs...)
13 +
    if method == :suave
14 +
        integrate = suave
15 +
    elseif method == :vegas
16 +
        integrate = vegas
17 +
    else
18 +
        ex = ErrorException("Integration method $method is not supported!")
19 +
        throw(ex)
20 +
    end
21 +
22 +
    units = unit(arg(x1)) * unit(x1)
23 +
24 +
    arg2(a) = ustrip(units, (x2 - x1) * arg((x2 - x1) * a + x1))::Float64
25 +
26 +
    function integrand(x, f)
27 +
        f[1] = arg2(x[1])
28 +
    end
29 +
30 +
    result, err = integrate(integrand, 1, 1; kwargs...)
31 +
32 +
    if units == Unitful.NoUnits
33 +
        return result[1], err[1]
34 +
    else
35 +
        return (result[1], err[1]) .* units
36 +
    end
37 +
end
38 +
39 +
#order, method, optional
40 +
41 +
function get_weights_fastgaussquadrature(order::Int, method::Symbol, optional::Real=NaN)
42 +
    if method == :gausslegendre
43 +
        xw = gausslegendre
44 +
        return xw(order)
45 +
46 +
    elseif method == :gausshermite
47 +
        xw = gausshermite
48 +
        return xw(order)
49 +
50 +
    elseif method == :gausslaguerre
51 +
        xw = gausslaguerre
52 +
        if isnan(optional)
53 +
            return xw(order)
54 +
        else
55 +
            α = optional
56 +
            return xw(order, α)
57 +
        end
58 +
59 +
    elseif method == :gausschebyshev
60 +
        xw = gausschebyshev
61 +
        if isnan(optional)
62 +
            return xw(order)
63 +
        else
64 +
            kind = optional
65 +
            return xw(order, kind)
66 +
        end
67 +
68 +
        # elseif method == :gaussjacobi
69 +
        #     xw = gaussjacobi
70 +
        #     return xw(order)
71 +
72 +
    elseif method == :gaussradau
73 +
        xw = gaussradau
74 +
        return xw(order)
75 +
76 +
    elseif method == :gausslobatto
77 +
        xw = gausslobatto
78 +
        return xw(order)
79 +
80 +
    else
81 +
        ex = ErrorException("Integration method $method is not supported!")
82 +
        throw(ex)
83 +
    end
84 +
85 +
86 +
end
87 +
88 +
const lru_xw_fastgaussquadrature = LRU{
89 +
    Tuple{Int,Symbol,Real},
90 +
    Tuple{Vector{Float64},Vector{Float64}},
91 +
}(maxsize=Int(1e5))
92 +
93 +
function get_weights_fastgaussquadrature_cached(order::Int, method::Symbol, optional::Real=NaN)
94 +
    get!(lru_xw_fastgaussquadrature, (order, method, optional)) do
95 +
        get_weights_fastgaussquadrature(order, method, optional)
96 +
    end
97 +
end
98 +
99 +
function _fastgaussquadrature_wrapper_1D(arg::Function, x1, x2; method::Symbol, order::Int, multithreading::Bool=true, verbose::Bool=false, kwargs...)
100 +
    if !(method in [:gausslegendre, :gausschebyshev, :gaussradau, :gausslobatto])
101 +
        ex = ErrorException("Integration method $method is not supported!")
102 +
        throw(ex)
103 +
    end
104 +
105 +
    if haskey(kwargs, :α)
106 +
        optional = kwargs[:α]
107 +
    elseif haskey(kwargs, :kind)
108 +
        optional = kwargs[:kind]
109 +
    else
110 +
        optional = NaN
111 +
    end
112 +
113 +
    x, w = get_weights_fastgaussquadrature_cached(order, method, optional)
114 +
115 +
    b = x2
116 +
    a = x1
117 +
118 +
    b_a_2 = (b - a) / 2
119 +
    a_plus_b_2 = (a + b) / 2
120 +
121 +
    arg_gauss(x) = arg(b_a_2 * x + a_plus_b_2)
122 +
123 +
    if verbose
124 +
        p = Progress(order)
125 +
    end
126 +
127 +
    if multithreading && order - 1 > nthreads() && nthreads() > 1
128 +
        tmp = arg_gauss(x[1])
129 +
        knots = zeros(order) * unit(tmp)
130 +
        knots[1] = tmp
131 +
        if verbose
132 +
            next!(p)
133 +
            @threads for i in 2:order
134 +
                knots[i] = arg_gauss(x[i])
135 +
                next!(p)
136 +
            end
137 +
        else
138 +
            @threads for i in 2:order
139 +
                knots[i] = arg_gauss(x[i])
140 +
            end
141 +
        end
142 +
    else
143 +
        knots = arg_gauss.(x)
144 +
    end
145 +
146 +
    integral = b_a_2 * dot(w, knots)
147 +
    return integral, NaN
148 +
149 +
end
150 +
151 +
function _fastgaussquadrature_wrapper_1D(arg::Function; method::Symbol, order::Int, multithreading::Bool=true, verbose::Bool=false, kwargs...)
152 +
    if !(method in [:gausshermite, :gausslaguerre])
153 +
        ex = ErrorException("Integration method $method is not supported!")
154 +
        throw(ex)
155 +
    end
156 +
157 +
    if haskey(kwargs, :α)
158 +
        optional = kwargs[:α]
159 +
    elseif haskey(kwargs, :kind)
160 +
        optional = kwargs[:kind]
161 +
    else
162 +
        optional = NaN
163 +
    end
164 +
165 +
    x, w = get_weights_fastgaussquadrature_cached(order, method, optional)
166 +
167 +
    if verbose
168 +
        p = Progress(order)
169 +
    end
170 +
171 +
    if multithreading && order - 1 > nthreads() && nthreads() > 1
172 +
        tmp = arg(x[1])
173 +
        knots = zeros(order) * unit(tmp)
174 +
        knots[1] = tmp
175 +
        if verbose
176 +
            next!(p)
177 +
            @threads for i in 2:order
178 +
                knots[i] = arg(x[i])
179 +
                next!(p)
180 +
            end
181 +
        else
182 +
            @threads for i in 2:order
183 +
                knots[i] = arg(x[i])
184 +
            end
185 +
        end
186 +
    else
187 +
        knots = arg.(x)
188 +
    end
189 +
190 +
    integral = dot(w, knots)
191 +
    return integral, NaN
192 +
193 +
end
194 +
195 +
@doc raw"""
196 +
    quad(arg::Function, x1, x2[; method = :quadgk, oder=1000, multithreading=true, kwargs...])
197 +
198 +
Performs the integral ``\int_{x1}^{x2}f(x)dx``
199 +
200 +
Available integration methods:
201 +
- `:suave`
202 +
- `:vegas`
203 +
- `:quadgk`
204 +
205 +
There are several fixed order quadrature methods available based on [`FastGaussQuadrature`](https://github.com/JuliaApproximation/FastGaussQuadrature.jl)
206 +
207 +
- `:gausslegendre` 
208 +
- `:gausshermite` 
209 +
- `:gausslaguerre` 
210 +
- `:gausschebyshev` 
211 +
- `:gaussradau` 
212 +
- `:gausslobatto`
213 +
214 +
If a specific integration routine is not needed, it is suggested to use `:quadgk` (the default option).
215 +
For highly oscillating integrals, it is advised to use `:gausslegendre`with a high order (~10000).
216 +
217 +
Please note that `:gausshermite` and `:gausslaguerre` are used to specific kind of integrals with infinite bounds.
218 +
219 +
For more informations on these integration techniques, please follow the official documentation [here](https://juliaapproximation.github.io/FastGaussQuadrature.jl/stable/gaussquadrature/).
220 +
221 +
See [QuadGK](https://github.com/JuliaMath/QuadGK.jl) and [Cuba.jl](https://giordano.github.io/Cuba.jl/stable/) for all the available keyword arguments.
222 +
# Examples
223 +
```jldoctest
224 +
using MultiQuad
225 +
226 +
func(x) = x^2*exp(-x)
227 +
quad(func, 0, 4) # equivalent to quad(func, 0, 4, method=:quadgk)
228 +
quad(func, 0, 4, method=:gausslegendre, order=1000)
229 +
230 +
# for certain kinds of integrals with infinite bounds, it may be possible to use a specific integration routine
231 +
func(x) = x^2*exp(-x)
232 +
quad(func, 0, Inf, method=:quadgk)[1]
233 +
# gausslaguerre computes the integral of f(x)*exp(-x) from 0 to infinity
234 +
quad(x -> x^4, method=:gausslaguerre, order=10000)[1] 
235 +
```
236 +
237 +
`quad` can handle units of measurement through the [`Unitful`](https://github.com/PainterQubits/Unitful.jl) package:
238 +
239 +
```jldoctest
240 +
using Unitful
241 +
242 +
func(x) = x^2
243 +
integral, error = quad(func, 1u"m", 5u"m")
244 +
```
245 +
"""
246 +
function quad(arg::Function, x1, x2; method::Symbol=:quadgk, order::Int=1000, multithreading::Bool=true, verbose::Bool=false, kwargs...)
247 +
248 +
    if method in [:suave, :vegas]
249 +
        return _cuba_wrapper_1D(arg, x1, x2; method=method)
250 +
    elseif method == :quadgk
251 +
        return _quadgk_wrapper_1D(arg, x1, x2; kwargs...)
252 +
    elseif method in [:gausslegendre, :gausschebyshev, :gaussradau, :gausslobatto]
253 +
        return _fastgaussquadrature_wrapper_1D(arg, x1, x2; method=method, order=order, multithreading=multithreading, verbose=verbose, kwargs...)
254 +
    else
255 +
        ex = ErrorException("Integration method $method is not supported!")
256 +
        throw(ex)
257 +
    end
258 +
end
259 +
260 +
@doc raw"""
261 +
    quad(arg::Function; method[, oder=1000, kwargs...])
262 +
"""
263 +
function quad(arg::Function; method::Symbol, order::Int=1000, multithreading::Bool=true, verbose::Bool=false, kwargs...)
264 +
265 +
    if method in [:gausshermite, :gausslaguerre]
266 +
        return _fastgaussquadrature_wrapper_1D(arg; method=method, order=order, multithreading=multithreading, verbose=verbose, kwargs...)
267 +
    else
268 +
        ex = ErrorException("Integration method $method is not supported!")
269 +
        throw(ex)
270 +
    end
271 +
end
272 +
273 +
function _cuba_wrapper_2D(arg::Function, x1, x2, y1::Function, y2::Function; method::Symbol, kwargs...)
274 +
    if method == :cuhre
275 +
        integrate = cuhre
276 +
    elseif method == :divonne
277 +
        integrate = divonne
278 +
    elseif method == :suave
279 +
        integrate = suave
280 +
    elseif method == :vegas
281 +
        integrate = vegas
282 +
    else
283 +
        ex = ErrorException("Integration method $method is not supported!")
284 +
        throw(ex)
285 +
    end
286 +
287 +
    units = unit(arg(y1(x1), x1)) * unit(x1) * unit(y1(x1))
288 +
289 +
    arg1(a, x) = (y2(x) - y1(x)) * arg((y2(x) - y1(x)) * a + y1(x), x)
290 +
291 +
292 +
    arg2(a, b) = ustrip(units, (x2 - x1) * arg1(a, (x2 - x1) * b + x1))::Float64
293 +
294 +
    function integrand2(x, f)
295 +
        f[1] = arg2(x[1], x[2])
296 +
    end
297 +
298 +
    result, err = integrate(integrand2, 2, 1; kwargs...)
299 +
300 +
    if units == Unitful.NoUnits
301 +
        return result[1], err[1]
302 +
    else
303 +
        return (result[1], err[1]) .* units
304 +
    end
305 +
end
306 +
307 +
function _cuba_wrapper_2D(arg::Function, x1, x2, y1, y2; method::Symbol, kwargs...)
308 +
    if method == :cuhre
309 +
        integrate = cuhre
310 +
    elseif method == :divonne
311 +
        integrate = divonne
312 +
    elseif method == :suave
313 +
        integrate = suave
314 +
    elseif method == :vegas
315 +
        integrate = vegas
316 +
    else
317 +
        ex = ErrorException("Integration method $method is not supported!")
318 +
        throw(ex)
319 +
    end
320 +
321 +
    units = unit(arg(y1, x1)) * unit(x1) * unit(y1)
322 +
323 +
    arg1(a, x) = (y2 - y1) * arg((y2 - y1) * a + y1, x)
324 +
325 +
    arg2(a, b) = ustrip(units, (x2 - x1) * arg1(a, (x2 - x1) * b + x1))::Float64
326 +
327 +
    function integrand2(x, f)
328 +
        f[1] = arg2(x[1], x[2])
329 +
    end
330 +
331 +
    result, err = integrate(integrand2, 2, 1; kwargs...)
332 +
333 +
    if units == Unitful.NoUnits
334 +
        return result[1], err[1]
335 +
    else
336 +
        return (result[1], err[1]) .* units
337 +
    end
338 +
end
339 +
340 +
341 +
function _hcubature_wrapper_2D(arg::Function, x1, x2, y1::Function, y2::Function; kwargs...)
342 +
    units = unit(arg(y1(x1), x1)) * unit(x1) * unit(y1(x1))
343 +
344 +
    arg1(a, x) = (y2(x) - y1(x)) * arg((y2(x) - y1(x)) * a + y1(x), x)
345 +
346 +
347 +
    arg2(a, b) = ustrip(units, (x2 - x1) * arg1(a, (x2 - x1) * b + x1))::Float64
348 +
349 +
350 +
    function integrand(arr)
351 +
        return arg2(arr[1], arr[2])
352 +
    end
353 +
354 +
    min_arr = [0, 0]
355 +
    max_arr = [1, 1]
356 +
    result, err = hcubature(integrand, min_arr, max_arr; kwargs...)
357 +
358 +
359 +
    if units == Unitful.NoUnits
360 +
        return result[1], err[1]
361 +
    else
362 +
        return (result[1], err[1]) .* units
363 +
    end
364 +
end
365 +
366 +
function _hcubature_wrapper_2D(arg::Function, x1, x2, y1, y2; kwargs...)
367 +
    units = unit(arg(y1, x1)) * unit(x1) * unit(y1)
368 +
369 +
    arg1(a, x) = (y2 - y1) * arg((y2 - y1) * a + y1, x)
370 +
371 +
    arg2(a, b) = ustrip(units, (x2 - x1) * arg1(a, (x2 - x1) * b + x1))::Float64
372 +
373 +
374 +
    function integrand(arr)
375 +
        return arg2(arr[1], arr[2])
376 +
    end
377 +
378 +
    min_arr = [0, 0]
379 +
    max_arr = [1, 1]
380 +
    result, err = hcubature(integrand, min_arr, max_arr; kwargs...)
381 +
382 +
383 +
    if units == Unitful.NoUnits
384 +
        return result[1], err[1]
385 +
    else
386 +
        return (result[1], err[1]) .* units
387 +
    end
388 +
end
389 +
390 +
@doc raw"""
391 +
    dblquad(arg::Function, x1, x2, y1::Function, y2::Function; method = :cuhre, kwargs...)
392 +
393 +
Performs the integral ``\int_{x1}^{x2}\int_{y1(x)}^{y2(x)}f(y,x)dydx``
394 +
395 +
Available integration methods:
396 +
- `:cuhre`
397 +
- `:divonne`
398 +
- `:suave`
399 +
- `:vegas`
400 +
- `:hcubature`
401 +
402 +
See [Cuba.jl](https://giordano.github.io/Cuba.jl/stable/) for all the available keyword arguments fro the `:cuhre`, `:divonne`, `:suave` and `:vegas` methods.
403 +
See [HCubature](https://github.com/stevengj/HCubature.jl) for all the available keywords for the `:hcubature` method.
404 +
405 +
# Examples
406 +
```jldoctest
407 +
func(y,x) = sin(x)*y^2
408 +
integral, error = dblquad(func, 1, 2, x->0, x->x^2, rtol=1e-9)
409 +
```
410 +
411 +
`dblquad` can handle units of measurement through the [`Unitful`](https://github.com/PainterQubits/Unitful.jl) package:
412 +
413 +
```jldoctest
414 +
using Unitful
415 +
416 +
func(y,x) = x*y^2
417 +
integral, error = dblquad(func, 1u"m", 2u"m", x->0u"m^2", x->x^2, rtol=1e-9)
418 +
```
419 +
"""
420 +
function dblquad(
421 +
    arg::Function,
422 +
    x1,
423 +
    x2,
424 +
    y1::Function,
425 +
    y2::Function;
426 +
    method=:hcubature,
427 +
    kwargs...
428 +
)
429 +
430 +
    if method in [:cuhre, :divonne, :suave, :vegas]
431 +
        return _cuba_wrapper_2D(arg, x1, x2, y1, y2; method=method)
432 +
    elseif method == :hcubature
433 +
        return _hcubature_wrapper_2D(arg, x1, x2, y1, y2; kwargs...)
434 +
    else
435 +
        ex = ErrorException("Integration method $method is not supported!")
436 +
        throw(ex)
437 +
    end
438 +
end
439 +
440 +
function dblquad(
441 +
    arg::Function,
442 +
    x1,
443 +
    x2,
444 +
    y1,
445 +
    y2;
446 +
    method=:hcubature,
447 +
    kwargs...
448 +
)
449 +
450 +
    if method in [:cuhre, :divonne, :suave, :vegas]
451 +
        return _cuba_wrapper_2D(arg, x1, x2, y1, y2; method=method)
452 +
    elseif method == :hcubature
453 +
        return _hcubature_wrapper_2D(arg, x1, x2, y1, y2; kwargs...)
454 +
    else
455 +
        ex = ErrorException("Integration method $method is not supported!")
456 +
        throw(ex)
457 +
    end
458 +
end
459 +
460 +
function _cuba_wrapper_3D(arg::Function, x1, x2, y1::Function, y2::Function, z1::Function, z2::Function; method::Symbol, kwargs...)
461 +
    if method == :cuhre
462 +
        integrate = cuhre
463 +
    elseif method == :divonne
464 +
        integrate = divonne
465 +
    elseif method == :suave
466 +
        integrate = suave
467 +
    elseif method == :vegas
468 +
        integrate = vegas
469 +
    else
470 +
        ex = ErrorException("Integration method $method is not supported!")
471 +
        throw(ex)
472 +
    end
473 +
474 +
    units = unit(arg(z1(x1, y1(x1)), y1(x1), x1)) * unit(x1) * unit(y1(x1)) *
475 +
            unit(z1(y1(x1), x1))
476 +
477 +
    arg0(a, y, x) =
478 +
        (z2(x, y) - z1(x, y)) * arg((z2(x, y) - z1(x, y)) * a + z1(x, y), y, x)
479 +
480 +
    arg1(a, b, x) = (y2(x) - y1(x)) * arg0(a, (y2(x) - y1(x)) * b + y1(x), x)
481 +
482 +
483 +
    arg2(a, b, c) =
484 +
        ustrip(units, (x2 - x1) * arg1(a, b, (x2 - x1) * c + x1))::Float64
485 +
486 +
487 +
    function integrand2(x, f)
488 +
        f[1] = arg2(x[1], x[2], x[3])
489 +
    end
490 +
491 +
    result, err = integrate(integrand2, 3, 1; kwargs...)
492 +
493 +
494 +
    if units == Unitful.NoUnits
495 +
        return result[1], err[1]
496 +
    else
497 +
        return (result[1], err[1]) .* units
498 +
    end
499 +
end
500 +
501 +
function _cuba_wrapper_3D(arg::Function, x1, x2, y1, y2, z1, z2; method::Symbol, kwargs...)
502 +
    if method == :cuhre
503 +
        integrate = cuhre
504 +
    elseif method == :divonne
505 +
        integrate = divonne
506 +
    elseif method == :suave
507 +
        integrate = suave
508 +
    elseif method == :vegas
509 +
        integrate = vegas
510 +
    else
511 +
        ex = ErrorException("Integration method $method is not supported!")
512 +
        throw(ex)
513 +
    end
514 +
515 +
    units = unit(arg(z1, y1, x1)) * unit(x1) * unit(y1) * unit(z1)
516 +
517 +
    arg0(a, y, x) = (z2 - z1) * arg((z2 - z1) * a + z1, y, x)
518 +
519 +
    arg1(a, b, x) = (y2 - y1) * arg0(a, (y2 - y1) * b + y1, x)
520 +
521 +
522 +
    arg2(a, b, c) =
523 +
        ustrip(units, (x2 - x1) * arg1(a, b, (x2 - x1) * c + x1))::Float64
524 +
525 +
526 +
    function integrand2(x, f)
527 +
        f[1] = arg2(x[1], x[2], x[3])
528 +
    end
529 +
530 +
    result, err = integrate(integrand2, 3, 1; kwargs...)
531 +
532 +
533 +
    if units == Unitful.NoUnits
534 +
        return result[1], err[1]
535 +
    else
536 +
        return (result[1], err[1]) .* units
537 +
    end
538 +
end
539 +
540 +
function _hcubature_wrapper_3D(arg::Function, x1, x2, y1::Function, y2::Function, z1::Function, z2::Function; kwargs...)
541 +
542 +
    units = unit(arg(z1(x1, y1(x1)), y1(x1), x1)) * unit(x1) * unit(y1(x1)) *
543 +
            unit(z1(y1(x1), x1))
544 +
545 +
    arg0(a, y, x) =
546 +
        (z2(x, y) - z1(x, y)) * arg((z2(x, y) - z1(x, y)) * a + z1(x, y), y, x)
547 +
548 +
    arg1(a, b, x) = (y2(x) - y1(x)) * arg0(a, (y2(x) - y1(x)) * b + y1(x), x)
549 +
550 +
551 +
    arg2(a, b, c) =
552 +
        ustrip(units, (x2 - x1) * arg1(a, b, (x2 - x1) * c + x1))::Float64
553 +
554 +
555 +
    function integrand(arr)
556 +
        return arg2(arr[1], arr[2], arr[3])
557 +
    end
558 +
559 +
    min_arr = [0, 0, 0]
560 +
    max_arr = [1, 1, 1]
561 +
    result, err = hcubature(integrand, min_arr, max_arr; kwargs...)
562 +
563 +
    if units == Unitful.NoUnits
564 +
        return result[1], err[1]
565 +
    else
566 +
        return (result[1], err[1]) .* units
567 +
    end
568 +
end
569 +
570 +
function _hcubature_wrapper_3D(arg::Function, x1, x2, y1, y2, z1, z2; kwargs...)
571 +
572 +
    units = unit(arg(z1, y1, x1)) * unit(x1) * unit(y1) * unit(z1)
573 +
574 +
    arg0(a, y, x) = (z2 - z1) * arg((z2 - z1) * a + z1, y, x)
575 +
576 +
    arg1(a, b, x) = (y2 - y1) * arg0(a, (y2 - y1) * b + y1, x)
577 +
578 +
579 +
    arg2(a, b, c) =
580 +
        ustrip(units, (x2 - x1) * arg1(a, b, (x2 - x1) * c + x1))::Float64
581 +
582 +
583 +
    function integrand(arr)
584 +
        return arg2(arr[1], arr[2], arr[3])
585 +
    end
586 +
587 +
    min_arr = [0, 0, 0]
588 +
    max_arr = [1, 1, 1]
589 +
    result, err = hcubature(integrand, min_arr, max_arr; kwargs...)
590 +
591 +
    if units == Unitful.NoUnits
592 +
        return result[1], err[1]
593 +
    else
594 +
        return (result[1], err[1]) .* units
595 +
    end
596 +
end
597 +
598 +
@doc raw"""
599 +
    tplquad(arg::Function, x1, x2, y1::Function, y2::Function, z1::Function, z2::Function; method = :cuhre, kwargs...)
600 +
601 +
Performs the integral ``\int_{x1}^{x2}\int_{y1(x)}^{y2(x)}\int_{z1(x,y)}^{z2(x,y)}f(z,y,x)dzdydx``
602 +
603 +
Available integration methods:
604 +
- `:cuhre`
605 +
- `:divonne`
606 +
- `:suave`
607 +
- `:vegas`
608 +
- `:hcubature`
609 +
610 +
See [Cuba.jl](https://giordano.github.io/Cuba.jl/stable/) for all the available keyword arguments fro the `:cuhre`, `:divonne`, `:suave` and `:vegas` methods.
611 +
See [HCubature](https://github.com/stevengj/HCubature.jl) for all the available keywords for the `:hcubature` method.
612 +
613 +
# Examples
614 +
```jldoctest
615 +
func(z,y,x) = sin(z)*y*x
616 +
integral, error = tplquad(func, 0, 4, x->x, x->x^2, (x,y)->2, (x,y)->3*x)
617 +
```
618 +
619 +
`tplquad` can handle units of measurement through the [`Unitful`](https://github.com/PainterQubits/Unitful.jl) package:
620 +
621 +
```jldoctest
622 +
using Unitful
623 +
624 +
func(z,y,x) = sin(z)*y*x
625 +
integral, error = tplquad(func, 0u"m", 4u"m", x->0u"m^2", x->x^2, (x,y)->0, (x,y)->3)
626 +
```
627 +
"""
628 +
function tplquad(
629 +
    arg::Function,
630 +
    x1,
631 +
    x2,
632 +
    y1::Function,
633 +
    y2::Function,
634 +
    z1::Function,
635 +
    z2::Function;
636 +
    method=:hcubature,
637 +
    kwargs...
638 +
)
639 +
640 +
    if method in [:cuhre, :divonne, :suave, :vegas]
641 +
        return _cuba_wrapper_3D(arg, x1, x2, y1, y2, z1, z2; method=method, kwargs...)
642 +
643 +
    elseif method == :hcubature
644 +
        return _hcubature_wrapper_3D(arg, x1, x2, y1, y2, z1, z2; kwargs...)
645 +
646 +
    else
647 +
        ex = ErrorException("Integration method $method is not supported!")
648 +
        throw(ex)
649 +
    end
650 +
end
651 +
652 +
function tplquad(
653 +
    arg::Function,
654 +
    x1,
655 +
    x2,
656 +
    y1,
657 +
    y2,
658 +
    z1,
659 +
    z2;
660 +
    method=:hcubature,
661 +
    kwargs...
662 +
)
663 +
664 +
    if method in [:cuhre, :divonne, :suave, :vegas]
665 +
        return _cuba_wrapper_3D(arg, x1, x2, y1, y2, z1, z2; method=method, kwargs...)
666 +
667 +
    elseif method == :hcubature
668 +
        return _hcubature_wrapper_3D(arg, x1, x2, y1, y2, z1, z2; kwargs...)
669 +
670 +
    else
671 +
        ex = ErrorException("Integration method $method is not supported!")
672 +
        throw(ex)
673 +
    end
674 +
end
675 +
676 +
677 +
end # module
Files Coverage
src/MultiQuad.jl 81.28%
Project Totals (1 files) 81.28%
Sunburst
The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.
Icicle
The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.
Grid
Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.
Loading