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 48
    if method == :suave
31 48
        integrate = suave
32 48
    elseif method == :vegas
33 48
        integrate = vegas
34 48
    elseif method == :quadgk
35 48
        integrate = quadgk
36
    else
37 36
        ex = ErrorException("Integration method $method is not supported!")
38 48
        throw(ex)
39
    end
40

41 48
    if method == :quadgk
42 48
        return quadgk(arg, x1, x2; kwargs...)
43
    end
44

45 36
    units = unit(arg(x1)) * unit(x1)
46

47 48
    arg2(a) = ustrip(units, (x2 - x1) * arg((x2 - x1) * a + x1))::Float64
48

49 48
    function integrand(x, f)
50 48
        f[1] = arg2(x[1])
51
    end
52

53 48
    result, err = integrate(integrand, 1, 1; kwargs...)
54

55 12
    if units == Unitful.NoUnits
56 48
        return result[1], err[1]
57
    else
58 48
        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 48
    if method == :cuhre
94 48
        integrate = cuhre
95 48
    elseif method == :divonne
96 48
        integrate = divonne
97 48
    elseif method == :suave
98 48
        integrate = suave
99 48
    elseif method == :vegas
100 48
        integrate = vegas
101
    else
102 36
        ex = ErrorException("Integration method $method is not supported!")
103 48
        throw(ex)
104
    end
105

106 48
    units = unit(arg(y1(x1), x1)) * unit(x1) * unit(y1(x1))
107

108 48
    arg1(a, x) = (y2(x) - y1(x)) * arg((y2(x) - y1(x)) * a + y1(x), x)
109

110

111 48
    arg2(a, b) = ustrip(units, (x2 - x1) * arg1(a, (x2 - x1) * b + x1))::Float64
112

113 48
    function integrand(x, f)
114 48
        f[1] = arg2(x[1], x[2])
115
    end
116

117 48
    result, err = integrate(integrand, 2, 1; kwargs...)
118

119 12
    if units == Unitful.NoUnits
120 48
        return result[1], err[1]
121
    else
122 48
        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 48
    if method == :cuhre
160 48
        integrate = cuhre
161 48
    elseif method == :divonne
162 48
        integrate = divonne
163 48
    elseif method == :suave
164 48
        integrate = suave
165 48
    elseif method == :vegas
166 48
        integrate = vegas
167
    else
168 36
        ex = ErrorException("Integration method $method is not supported!")
169 48
        throw(ex)
170
    end
171

172 48
    units = unit(arg(z1(x1, y1(x1)), y1(x1), x1)) * unit(x1) * unit(y1(x1)) *
173
            unit(z1(y1(x1), x1))
174

175 48
    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 48
    arg1(a, b, x) = (y2(x) - y1(x)) * arg0(a, (y2(x) - y1(x)) * b + y1(x), x)
179

180

181 48
    arg2(a, b, c) =
182
        ustrip(units, (x2 - x1) * arg1(a, b, (x2 - x1) * c + x1))::Float64
183

184 48
    function integrand(x, f)
185 48
        f[1] = arg2(x[1], x[2], x[3])
186
    end
187

188 48
    result, err = integrate(integrand, 3, 1; kwargs...)
189

190 12
    if units == Unitful.NoUnits
191 48
        return result[1], err[1]
192
    else
193 48
        return (result[1], err[1]) .* units
194
    end
195
end
196

197

198
end # module

Read our documentation on viewing source code .

Loading