1
module MultiQuad
2

3
using QuadGK, Cuba, HCubature, 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

12
Available integration methods:
13
- `:suave`
14
- `:vegas`
15
- `:quadgk`
16

17
See [QuadGK](https://github.com/JuliaMath/QuadGK.jl) and [Cuba.jl](https://giordano.github.io/Cuba.jl/stable/) for all the available keyword arguments.
18

19
# Examples
20
```jldoctest
21
func(x) = x^2*exp(-x)
22
integral, error = quad(func, 0, 4)
23
```
24

25
`quad` can handle units of measurement through the [`Unitful`](https://github.com/PainterQubits/Unitful.jl) package:
26

27
```jldoctest
28
using Unitful
29

30
func(x) = x^2
31
integral, error = quad(func, 1u"m", 5u"m")
32
```
33
"""
34
function quad(arg::Function, x1, x2; method = :quadgk, kwargs...)
35

36 30
    if method == :suave
37 30
        integrate = suave
38 30
    elseif method == :vegas
39 30
        integrate = vegas
40 30
    elseif method == :quadgk
41 30
        integrate = quadgk
42
    else
43 24
        ex = ErrorException("Integration method $method is not supported!")
44 30
        throw(ex)
45
    end
46

47 30
    if method == :quadgk
48 30
        return quadgk(arg, x1, x2; kwargs...)
49
    end
50

51 24
    units = unit(arg(x1)) * unit(x1)
52

53 30
    arg2(a) = ustrip(units, (x2 - x1) * arg((x2 - x1) * a + x1))::Float64
54

55 30
    function integrand(x, f)
56 30
        f[1] = arg2(x[1])
57
    end
58

59 30
    result, err = integrate(integrand, 1, 1; kwargs...)
60

61 12
    if units == Unitful.NoUnits
62 30
        return result[1], err[1]
63
    else
64 30
        return (result[1], err[1]) .* units
65
    end
66
end
67

68
@doc raw"""
69
    dblquad(arg::Function, x1, x2, y1::Function, y2::Function; method = :cuhre, kwargs...)
70

71
Performs the integral ``\int_{x1}^{x2}\int_{y1(x)}^{y2(x)}f(y,x)dydx``
72

73
Available integration methods:
74
- `:cuhre`
75
- `:divonne`
76
- `:suave`
77
- `:vegas`
78
- `:hcubature`
79

80
See [Cuba.jl](https://giordano.github.io/Cuba.jl/stable/) for all the available keyword arguments fro the `:cuhre`, `:divonne`, `:suave` and `:vegas` methods.
81
See [HCubature](https://github.com/stevengj/HCubature.jl) for all the available keywords for the `:hcubature` method.
82

83
# Examples
84
```jldoctest
85
func(y,x) = sin(x)*y^2
86
integral, error = dblquad(func, 1, 2, x->0, x->x^2, rtol=1e-9)
87
```
88

89
`dblquad` can handle units of measurement through the [`Unitful`](https://github.com/PainterQubits/Unitful.jl) package:
90

91
```jldoctest
92
using Unitful
93

94
func(y,x) = x*y^2
95
integral, error = dblquad(func, 1u"m", 2u"m", x->0u"m^2", x->x^2, rtol=1e-9)
96
```
97
"""
98
function dblquad(
99
    arg::Function,
100
    x1,
101
    x2,
102
    y1::Function,
103
    y2::Function;
104
    method = :cuhre,
105
    kwargs...,
106
)
107

108 30
    if method == :cuhre
109 30
        integrate = cuhre
110 30
    elseif method == :divonne
111 30
        integrate = divonne
112 30
    elseif method == :suave
113 30
        integrate = suave
114 30
    elseif method == :vegas
115 30
        integrate = vegas
116 30
    elseif method == :hcubature
117 30
        integrate = hcubature
118
    else
119 24
        ex = ErrorException("Integration method $method is not supported!")
120 30
        throw(ex)
121
    end
122

123 30
    units = unit(arg(y1(x1), x1)) * unit(x1) * unit(y1(x1))
124

125 30
    arg1(a, x) = (y2(x) - y1(x)) * arg((y2(x) - y1(x)) * a + y1(x), x)
126

127

128 30
    arg2(a, b) = ustrip(units, (x2 - x1) * arg1(a, (x2 - x1) * b + x1))::Float64
129

130 30
    if method == :hcubature
131 30
        function integrand(arr)
132 30
            return arg2(arr[1], arr[2])
133
        end
134

135 30
        min_arr = [0, 0]
136 30
        max_arr = [1, 1]
137 30
        result, err = integrate(integrand, min_arr, max_arr; kwargs...)
138
    else
139 30
        function integrand2(x, f)
140 30
            f[1] = arg2(x[1], x[2])
141
        end
142

143 30
        result, err = integrate(integrand2, 2, 1; kwargs...)
144
    end
145

146 30
    if units == Unitful.NoUnits
147 30
        return result[1], err[1]
148
    else
149 30
        return (result[1], err[1]) .* units
150
    end
151
end
152

153
@doc raw"""
154
    tplquad(arg::Function, x1, x2, y1::Function, y2::Function, z1::Function, z2::Function; method = :cuhre, kwargs...)
155

156
Performs the integral ``\int_{x1}^{x2}\int_{y1(x)}^{y2(x)}\int_{z1(x,y)}^{z2(x,y)}f(z,y,x)dzdydx``
157

158
Available integration methods:
159
- `:cuhre`
160
- `:divonne`
161
- `:suave`
162
- `:vegas`
163
- `:hcubature`
164

165
See [Cuba.jl](https://giordano.github.io/Cuba.jl/stable/) for all the available keyword arguments fro the `:cuhre`, `:divonne`, `:suave` and `:vegas` methods.
166
See [HCubature](https://github.com/stevengj/HCubature.jl) for all the available keywords for the `:hcubature` method.
167

168
# Examples
169
```jldoctest
170
func(z,y,x) = sin(z)*y*x
171
integral, error = tplquad(func, 0, 4, x->x, x->x^2, (x,y)->2, (x,y)->3*x)
172
```
173

174
`tplquad` can handle units of measurement through the [`Unitful`](https://github.com/PainterQubits/Unitful.jl) package:
175

176
```jldoctest
177
using Unitful
178

179
func(z,y,x) = sin(z)*y*x
180
integral, error = tplquad(func, 0u"m", 4u"m", x->0u"m^2", x->x^2, (x,y)->0, (x,y)->3)
181
```
182
"""
183
function tplquad(
184
    arg::Function,
185
    x1,
186
    x2,
187
    y1::Function,
188
    y2::Function,
189
    z1::Function,
190
    z2::Function;
191
    method = :cuhre,
192
    kwargs...,
193
)
194

195 30
    if method == :cuhre
196 30
        integrate = cuhre
197 30
    elseif method == :divonne
198 30
        integrate = divonne
199 30
    elseif method == :suave
200 30
        integrate = suave
201 30
    elseif method == :vegas
202 30
        integrate = vegas
203 30
    elseif method == :hcubature
204 30
        integrate = hcubature
205
    else
206 24
        ex = ErrorException("Integration method $method is not supported!")
207 30
        throw(ex)
208
    end
209

210 30
    units = unit(arg(z1(x1, y1(x1)), y1(x1), x1)) * unit(x1) * unit(y1(x1)) *
211
            unit(z1(y1(x1), x1))
212

213 30
    arg0(a, y, x) =
214
        (z2(x, y) - z1(x, y)) * arg((z2(x, y) - z1(x, y)) * a + z1(x, y), y, x)
215

216 30
    arg1(a, b, x) = (y2(x) - y1(x)) * arg0(a, (y2(x) - y1(x)) * b + y1(x), x)
217

218

219 30
    arg2(a, b, c) =
220
        ustrip(units, (x2 - x1) * arg1(a, b, (x2 - x1) * c + x1))::Float64
221

222 30
    if method == :hcubature
223 30
        function integrand(arr)
224 30
            return arg2(arr[1], arr[2], arr[3])
225
        end
226

227 30
        min_arr = [0, 0, 0]
228 30
        max_arr = [1, 1, 1]
229 30
        result, err = integrate(integrand, min_arr, max_arr; kwargs...)
230
    else
231 30
        function integrand2(x, f)
232 30
            f[1] = arg2(x[1], x[2], x[3])
233
        end
234

235 30
        result, err = integrate(integrand2, 3, 1; kwargs...)
236
    end
237

238 30
    if units == Unitful.NoUnits
239 30
        return result[1], err[1]
240
    else
241 30
        return (result[1], err[1]) .* units
242
    end
243
end
244

245

246
end # module

Read our documentation on viewing source code .

Loading