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 .

Loading