Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F63292412
IntegralSolver2D.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sun, May 19, 02:32
Size
6 KB
Mime Type
text/x-c
Expires
Tue, May 21, 02:32 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
17756054
Attached To
R7917 Guggisberg-Liardon - Numerical integration
IntegralSolver2D.cpp
View Options
#include "IntegralSolver2D.hpp"
#include <iostream>
// Constructor and destructor
IntegralSolver2D
::
IntegralSolver2D
(
FunctionRntoR
function_
,
DomainCartesian
domain_
,
Method
method_
)
:
IntegrationSolver
<
FunctionRntoR
,
DomainCartesian
>
(
function_
,
domain_
,
method_
){}
IntegralSolver2D
::~
IntegralSolver2D
()
{}
// Get result method
double
IntegralSolver2D
::
Compute
()
{
// Method that will be used
type
=
method
.
name
;
// Extracting the boundaries of the domain and subdivision
a
=
domain
.
boundaries
[
0
][
0
];
b
=
domain
.
boundaries
[
0
][
1
];
c
=
domain
.
boundaries
[
1
][
0
];
d
=
domain
.
boundaries
[
1
][
1
];
n
=
domain
.
subdivision
[
0
];
// Number of sub-intervals in the interval [a,b] of X
m
=
domain
.
subdivision
[
1
];
// Number of sub-intervals in the interval [c,d] of Y
// Computing of the result the method given by type
if
(
type
==
"Midpoint"
)
{
return
this
->
IntegrationMidpoint
();
}
if
(
type
==
"Trapz"
)
{
return
this
->
IntegrationTrapeze
();
}
if
(
type
==
"Simpson"
)
{
return
this
->
IntegrationSimpson
();
}
else
{
OutOfRangeException
error
(
"This integration method doesn't exist."
);
error
.
PrintDebug
();
exit
(
0
);
}
}
/* Midpoint method
*
* For each cell, the function value in the middle is computed and multiplied by the area of the cell.
*
*/
double
IntegralSolver2D
::
IntegrationMidpoint
()
{
double
result
=
0
;
double
sizeX
=
(
b
-
a
)
/
n
;
// Size of the cell on the X axis
double
sizeY
=
(
d
-
c
)
/
m
;
// Size of the cell on the Y axis
// Iteration on each cell,
for
(
int
i
=
0
;
i
<
n
;
i
++
){
for
(
int
j
=
0
;
j
<
m
;
j
++
){
// Position where to evaluate the function
vector
<
double
>
arg
;
arg
.
push_back
(
i
*
sizeX
+
sizeX
/
2
);
arg
.
push_back
(
j
*
sizeY
+
sizeY
/
2
);
result
+=
sizeX
*
sizeY
*
function
.
GetResult
(
arg
);
}
}
return
result
;
}
/* Trapeze method
*
* For each cell, takes the average function values from the corners and multiply it by the area of the cell
*
*/
double
IntegralSolver2D
::
IntegrationTrapeze
()
{
double
result
=
0
;
double
sizeX
=
(
b
-
a
)
/
n
;
// Size of the cell on the X axis
double
sizeY
=
(
d
-
c
)
/
m
;
// Size of the cell on the Y axis
// Iteration on each cell
for
(
int
i
=
0
;
i
<
n
;
i
++
){
for
(
int
j
=
0
;
j
<
m
;
j
++
)
{
// Positions of the corner where to evaluate the function
vector
<
double
>
corner1
;
corner1
.
push_back
(
i
*
sizeX
);
corner1
.
push_back
(
j
*
sizeY
);
vector
<
double
>
corner2
;
corner2
.
push_back
((
i
+
1
)
*
sizeX
);
corner2
.
push_back
(
j
*
sizeY
);
vector
<
double
>
corner3
;
corner3
.
push_back
(
i
*
sizeX
);
corner3
.
push_back
((
j
+
1
)
*
sizeY
);
vector
<
double
>
corner4
;
corner4
.
push_back
((
i
+
1
)
*
sizeX
);
corner4
.
push_back
((
j
+
1
)
*
sizeY
);
result
+=
sizeX
*
sizeY
*
(
function
.
GetResult
(
corner1
)
+
function
.
GetResult
(
corner2
)
+
function
.
GetResult
(
corner3
)
+
function
.
GetResult
(
corner4
))
/
4
;
}
}
return
result
;
}
/* Simpson method
*
* Algorithm taken from : http://mathfaculty.fullerton.edu/mathews/n2003/SimpsonsRule2DMod.html
*
* Uses a quadratic approximation of the function for each cell and compute its integral on the cell domain
*
*/
double
IntegralSolver2D
::
IntegrationSimpson
()
{
double
result
=
0
;
double
h
=
(
b
-
a
)
/
(
2
*
n
);
double
k
=
(
d
-
c
)
/
(
2
*
m
);
// Vector that will contains some of the position where to evaluate the function
vector
<
double
>
a_const
;
a_const
.
push_back
(
a
);
a_const
.
push_back
(
c
);
vector
<
double
>
c_const
;
c_const
.
push_back
(
a
);
c_const
.
push_back
(
d
);
vector
<
double
>
b_const
;
b_const
.
push_back
(
b
);
b_const
.
push_back
(
c
);
vector
<
double
>
d_const
;
d_const
.
push_back
(
a
);
d_const
.
push_back
(
d
);
result
+=
function
.
GetResult
(
a_const
)
+
function
.
GetResult
(
c_const
)
+
function
.
GetResult
(
b_const
)
+
function
.
GetResult
(
d_const
);
for
(
int
j
=
1
;
j
<=
n
;
j
++
){
a_const
[
1
]
=
c
+
(
2
*
j
-
1
)
*
k
;
b_const
[
1
]
=
c
+
(
2
*
j
-
1
)
*
k
;
result
+=
4
*
function
.
GetResult
(
a_const
)
+
4
*
function
.
GetResult
(
b_const
);
}
for
(
int
j
=
1
;
j
<=
n
-
1
;
j
++
){
a_const
[
1
]
=
c
+
(
2
*
j
)
*
k
;
b_const
[
1
]
=
c
+
(
2
*
j
)
*
k
;
result
+=
2
*
function
.
GetResult
(
a_const
)
+
2
*
function
.
GetResult
(
b_const
);
}
for
(
int
i
=
1
;
i
<=
m
;
i
++
){
c_const
[
0
]
=
a
+
(
2
*
i
-
1
)
*
k
;
d_const
[
0
]
=
a
+
(
2
*
i
-
1
)
*
k
;
result
+=
4
*
function
.
GetResult
(
c_const
)
+
4
*
function
.
GetResult
(
d_const
);
}
for
(
int
i
=
1
;
i
<=
m
-
1
;
i
++
){
c_const
[
0
]
=
a
+
(
2
*
i
)
*
k
;
d_const
[
0
]
=
a
+
(
2
*
i
)
*
k
;
result
+=
2
*
function
.
GetResult
(
c_const
)
+
2
*
function
.
GetResult
(
d_const
);
}
for
(
int
i
=
1
;
i
<=
n
;
i
++
){
for
(
int
j
=
1
;
j
<=
m
;
j
++
){
// Position where to evaluate the function
vector
<
double
>
arg
;
arg
.
push_back
(
a
+
(
2
*
i
-
1
)
*
h
);
arg
.
push_back
(
c
+
(
2
*
j
-
1
)
*
k
);
result
+=
16
*
function
.
GetResult
(
arg
);
}
}
for
(
int
i
=
1
;
i
<=
n
-
1
;
i
++
){
for
(
int
j
=
1
;
j
<=
m
;
j
++
){
// Position where to evaluate the function
vector
<
double
>
arg
;
arg
.
push_back
(
a
+
(
2
*
i
-
1
)
*
h
);
arg
.
push_back
(
c
+
(
2
*
j
)
*
k
);
result
+=
8
*
function
.
GetResult
(
arg
);
}
}
for
(
int
i
=
1
;
i
<=
n
;
i
++
){
for
(
int
j
=
1
;
j
<=
m
-
1
;
j
++
){
// Position where to evaluate the function
vector
<
double
>
arg
;
arg
.
push_back
(
a
+
(
2
*
i
)
*
h
);
arg
.
push_back
(
c
+
(
2
*
j
-
1
)
*
k
);
result
+=
8
*
function
.
GetResult
(
arg
);
}
}
for
(
int
i
=
1
;
i
<=
n
-
1
;
i
++
){
for
(
int
j
=
1
;
j
<=
m
-
1
;
j
++
){
// Position where to evaluate the function
vector
<
double
>
arg
;
arg
.
push_back
(
a
+
(
2
*
i
)
*
h
);
arg
.
push_back
(
c
+
(
2
*
j
)
*
k
);
result
+=
4
*
function
.
GetResult
(
arg
);
}
}
result
=
result
*
h
*
k
/
9
;
return
result
;
}
Event Timeline
Log In to Comment