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 2
    if method == :suave
37 2
        integrate = suave
38 2
    elseif method == :vegas
39 2
        integrate = vegas
40 2
    elseif method == :quadgk
41 2
        integrate = quadgk
42
    else
43 1
        ex = ErrorException("Integration method $method is not supported!")
44 2
        throw(ex)
45
    end
46

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

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

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

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

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

61 1
    if units == Unitful.NoUnits
62 2
        return result[1], err[1]
63
    else
64 2
        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 2
    if method == :cuhre
109 2
        integrate = cuhre
110 2
    elseif method == :divonne
111 2
        integrate = divonne
112 2
    elseif method == :suave
113 2
        integrate = suave
114 2
    elseif method == :vegas
115 2
        integrate = vegas
116 2
    elseif method == :hcubature
117 2
        integrate = hcubature
118
    else
119 1
        ex = ErrorException("Integration method $method is not supported!")
120 2
        throw(ex)
121
    end
122

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

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

127

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

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

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

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

146 2
    if units == Unitful.NoUnits
147 2
        return result[1], err[1]
148
    else
149 2
        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 2
    if method == :cuhre
196 2
        integrate = cuhre
197 2
    elseif method == :divonne
198 2
        integrate = divonne
199 2
    elseif method == :suave
200 2
        integrate = suave
201 2
    elseif method == :vegas
202 2
        integrate = vegas
203 2
    elseif method == :hcubature
204 2
        integrate = hcubature
205
    else
206 1
        ex = ErrorException("Integration method $method is not supported!")
207 2
        throw(ex)
208
    end
209

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

213 2
    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 2
    arg1(a, b, x) = (y2(x) - y1(x)) * arg0(a, (y2(x) - y1(x)) * b + y1(x), x)
217

218

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

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

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

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

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

245

246
end # module

Read our documentation on viewing source code .

Loading