|
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
|
17
|
if method == :suave
|
|
31
|
17
|
integrate = suave
|
|
32
|
17
|
elseif method == :vegas
|
|
33
|
17
|
integrate = vegas
|
|
34
|
17
|
elseif method == :quadgk
|
|
35
|
17
|
integrate = quadgk
|
|
36
|
|
else
|
|
37
|
13
|
ex = ErrorException("Integration method $method is not supported!")
|
|
38
|
17
|
throw(ex)
|
|
39
|
|
end
|
|
40
|
|
|
|
41
|
17
|
if method == :quadgk
|
|
42
|
17
|
return quadgk(arg, x1, x2; kwargs...)
|
|
43
|
|
end
|
|
44
|
|
|
|
45
|
13
|
units = unit(arg(x1)) * unit(x1)
|
|
46
|
|
|
|
47
|
17
|
arg2(a) = ustrip(units, (x2 - x1) * arg((x2 - x1) * a + x1))::Float64
|
|
48
|
|
|
|
49
|
17
|
function integrand(x, f)
|
|
50
|
17
|
f[1] = arg2(x[1])
|
|
51
|
|
end
|
|
52
|
|
|
|
53
|
17
|
result, err = integrate(integrand, 1, 1; kwargs...)
|
|
54
|
|
|
|
55
|
4
|
if units == Unitful.NoUnits
|
|
56
|
17
|
return result[1], err[1]
|
|
57
|
|
else
|
|
58
|
17
|
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
|
17
|
if method == :cuhre
|
|
94
|
17
|
integrate = cuhre
|
|
95
|
17
|
elseif method == :divonne
|
|
96
|
17
|
integrate = divonne
|
|
97
|
17
|
elseif method == :suave
|
|
98
|
17
|
integrate = suave
|
|
99
|
17
|
elseif method == :vegas
|
|
100
|
17
|
integrate = vegas
|
|
101
|
|
else
|
|
102
|
13
|
ex = ErrorException("Integration method $method is not supported!")
|
|
103
|
17
|
throw(ex)
|
|
104
|
|
end
|
|
105
|
|
|
|
106
|
17
|
units = unit(arg(y1(x1), x1)) * unit(x1) * unit(y1(x1))
|
|
107
|
|
|
|
108
|
17
|
arg1(a, x) = (y2(x) - y1(x)) * arg((y2(x) - y1(x)) * a + y1(x), x)
|
|
109
|
|
|
|
110
|
|
|
|
111
|
17
|
arg2(a, b) = ustrip(units, (x2 - x1) * arg1(a, (x2 - x1) * b + x1))::Float64
|
|
112
|
|
|
|
113
|
17
|
function integrand(x, f)
|
|
114
|
17
|
f[1] = arg2(x[1], x[2])
|
|
115
|
|
end
|
|
116
|
|
|
|
117
|
17
|
result, err = integrate(integrand, 2, 1; kwargs...)
|
|
118
|
|
|
|
119
|
4
|
if units == Unitful.NoUnits
|
|
120
|
17
|
return result[1], err[1]
|
|
121
|
|
else
|
|
122
|
17
|
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
|
17
|
if method == :cuhre
|
|
160
|
17
|
integrate = cuhre
|
|
161
|
17
|
elseif method == :divonne
|
|
162
|
17
|
integrate = divonne
|
|
163
|
17
|
elseif method == :suave
|
|
164
|
17
|
integrate = suave
|
|
165
|
17
|
elseif method == :vegas
|
|
166
|
17
|
integrate = vegas
|
|
167
|
|
else
|
|
168
|
13
|
ex = ErrorException("Integration method $method is not supported!")
|
|
169
|
17
|
throw(ex)
|
|
170
|
|
end
|
|
171
|
|
|
|
172
|
17
|
units = unit(arg(z1(x1, y1(x1)), y1(x1), x1)) * unit(x1) * unit(y1(x1)) *
|
|
173
|
|
unit(z1(y1(x1), x1))
|
|
174
|
|
|
|
175
|
17
|
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
|
17
|
arg1(a, b, x) = (y2(x) - y1(x)) * arg0(a, (y2(x) - y1(x)) * b + y1(x), x)
|
|
179
|
|
|
|
180
|
|
|
|
181
|
17
|
arg2(a, b, c) =
|
|
182
|
|
ustrip(units, (x2 - x1) * arg1(a, b, (x2 - x1) * c + x1))::Float64
|
|
183
|
|
|
|
184
|
17
|
function integrand(x, f)
|
|
185
|
17
|
f[1] = arg2(x[1], x[2], x[3])
|
|
186
|
|
end
|
|
187
|
|
|
|
188
|
17
|
result, err = integrate(integrand, 3, 1; kwargs...)
|
|
189
|
|
|
|
190
|
4
|
if units == Unitful.NoUnits
|
|
191
|
17
|
return result[1], err[1]
|
|
192
|
|
else
|
|
193
|
17
|
return (result[1], err[1]) .* units
|
|
194
|
|
end
|
|
195
|
|
end
|
|
196
|
|
|
|
197
|
|
|
|
198
|
|
end # module
|
KristofferC