-
Notifications
You must be signed in to change notification settings - Fork 2
/
Integration.jl
82 lines (75 loc) · 1.9 KB
/
Integration.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
"""
CalcSingleIntegral(f::Function, a::Real, b::Real, n::UInt64)
Compute integral of f function from a to b (Composite Simpson method).
b
⌠
⌡ f(x) ⋅ dx
a
"""
function CalcSingleIntegral(f::Function, a::Real, b::Real, n::UInt64)
xi = 0.0
h = (b - a) / n
xi0 = f(a) + f(b)
xi1 = 0.0
xi2 = 0.0
for i = 1:(n - 1)
x = a + i * h;
if mod(i, 2) == 0
xi2 = xi2 + f(x)
else
xi1 = xi1 + f(x)
end
end
xi = h * (xi0 + 2 * xi2 + 4 * xi1) / 3;
println("Approximation to integral: $xi")
return xi
end
"""
CalcDoubleIntegral(f::Function, a::Real, b::Real, c::Real, d::Real, n::UInt64)
Compute double integral of f(x,y) function over domain D=[a, b] x [c(a), d(b)]
(Double Simpson's method).
b
.- d(x)
| .-
| |
| -' f(x, y) . dy . dx
-' c(x)
a
"""
function CalcDoubleIntegral(f::Function, a::Real, b::Real, c::Function, d::Function, m::UInt64, n::UInt64)
@assert rem(m, 2) == 0
@assert rem(n, 2) == 0
hx = (b - a) / m
t = 0.0
for i = 0:m
x = a + (i * hx)
cx = c(x)
dx = d(x)
hy = (dx - cx) / n;
s = 0.0
for j = 0:n
y = cx + (j * hy)
fxy = f(x, y)
w = 4.0
if rem(j, 2) == 0
w = 2.0
end
if (j == 0) | (j == n)
w = 1.0
end
s = s + (w * fxy);
end
s = s * hy / 3.0
w = 4.0
if rem(i, 2) == 0
w = 2.0
end
if (i == 0) | (i == m)
w = 1.0
end
t = t + (w * s)
end
t = t * hx / 3.0
println("Approximation to integral: $t")
return t
end