Homec4science
Diffusion Guggisberg-Liardon - Numerical integration (master) Browse

R7917/IntegrationSolver2DCartesian.cpp
0f8ddfe1ee20master

Last Change
Raw File

IntegrationSolver2DCartesian.cpp
View Options

#include "IntegrationSolver2DCartesian.hpp"
#include <iostream>
// Constructor and destructor
IntegrationSolver2DCartesian::IntegrationSolver2DCartesian(FunctionRntoR function_,DomainCartesian domain_, Method method_)
: IntegrationSolver<FunctionRntoR, DomainCartesian, double>(function_, domain_, method_){}
IntegrationSolver2DCartesian::~IntegrationSolver2DCartesian() {}
// Compute method
double IntegrationSolver2DCartesian::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 IntegrationSolver2DCartesian::IntegrationMidpoint() {
double result = 0;
double sizeX = (b - a) / m; // Size of the cell on the X axis
double sizeY = (d - c) / n; // 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 IntegrationSolver2DCartesian::IntegrationTrapeze() {
double result = 0;
double sizeX = (b - a) / m; // Size of the cell on the X axis
double sizeY = (d - c) / n; // 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 IntegrationSolver2DCartesian::IntegrationSimpson() {
double result = 0;
double h = (b - a) / (2 * m);
double k = (d - c) / (2 * n);
// 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;
}
c4science · Help