Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <iostream>
- #include <cmath>
- #include <algorithm>
- #include <vector>
- #include <map>
- using namespace std;
- double x_left, x_right, y_bottom, y_top;
- double h_x, h_y;
- int M, N;
- vector<double> w;
- // exact solution u(x, y)
- double u(double x, double y) {
- return 2.0 / (1 + pow(x, 2) + pow(y, 2));
- }
- double k(double x, double y) {
- return 1 + pow(x + y, 2);
- }
- // potential function
- double q(double x, double y) {
- return 1;
- }
- // psi_R(x, y) = k(x_right, y) * du/dx(x_right, y)
- double psi_R(double y) {
- return (-12) * (pow((y + 3), 2) + 1) / pow((pow(y, 2) + 10), 2);
- }
- // psi_L(x, y) = -k(x_left, y) * du/dx(x_left, y)
- double psi_L(double y) {
- return (-8) * (pow((y - 2), 2) + 1) / pow((pow(y, 2) + 5), 2);
- }
- // psi_T(x, y) = k(x, y_top) * du/dy(x, y_top)
- double psi_T(double x) {
- return (-16) * (pow((x + 4), 2) + 1) / pow((pow(x, 2) + 17), 2);
- }
- // psi_B(x, y) = -k(x, y_top) * du/dy(x, y_top)
- double psi_B(double x) {
- return (-4) * (pow((x - 1), 2) + 1) / pow((pow(x, 2) + 2), 2);
- }
- // right-part function of Poisson equation
- double F(double x, double y) {
- return 2 * (pow(x,4) + pow(y,4) + 2 * (pow(x,2) + 3) * pow(y,2) + 6 * pow(x,2) + 16*x*y + 5)
- / pow((1 + pow(x, 2) + pow(y, 2)), 3);
- }
- double Get(int i, int j){
- return w[j*(M + 1) + i];
- }
- void CalculateExactSolution(vector<double> &exact) {
- for (int i = 0; i < M + 1; ++i)
- for (int j = 0; j < N + 1; ++j)
- exact[j * (M + 1) + i] = u(x_left + i * h_x, y_bottom + j * h_y);
- //cout << "Exact solution of Poisson equation:" << endl;
- //PrintVector(exact);
- //cout << endl;
- }
- int IsNodeInternalCornerOrSide(int current_node_global_offset_x, int current_node_global_offset_y){
- //corners
- //left bottom corner
- if (current_node_global_offset_x == 0 && current_node_global_offset_y == 0){
- return 2;
- }
- //left top corner
- if (current_node_global_offset_x == 0 && current_node_global_offset_y == N){
- return 4;
- }
- //right bottom corner
- if (current_node_global_offset_x == M && current_node_global_offset_y == 0){
- return 6;
- }
- //right top corner
- if (current_node_global_offset_x == M && current_node_global_offset_y == N){
- return 8;
- }
- //sides
- //left side
- if (current_node_global_offset_y >= 1 && current_node_global_offset_y <= N - 1 &&
- current_node_global_offset_x == 0){
- return 1;
- }
- //right side
- if (current_node_global_offset_y >= 1 && current_node_global_offset_y <= N - 1 &&
- current_node_global_offset_x == M){
- return 3;
- }
- //bottom side
- if (current_node_global_offset_x >= 1 && current_node_global_offset_x <= M - 1 &&
- current_node_global_offset_y == 0){
- return 5;
- }
- //top side
- if ((current_node_global_offset_x >= 1 && current_node_global_offset_x <= M - 1 &&
- current_node_global_offset_y == N)){
- return 7;
- }
- //internal
- if ((current_node_global_offset_x >= 1 && current_node_global_offset_x <= M - 1) &&
- (current_node_global_offset_y >= 1 && current_node_global_offset_y <= N - 1)){
- return 0;
- }
- return -1;
- }
- double ComputeMagicInnerProductA_iw (int current_node_global_offset_x, int current_node_global_offset_y){
- int glob_x = current_node_global_offset_x;
- int glob_y = current_node_global_offset_y;
- double result = 0.0;
- map <string,bool> neighbours = {
- {"left", true},
- {"right", true},
- {"bottom", true},
- {"top", true}
- };
- double left_neighbour = 0.0, right_neighbour = 0.0, bottom_neighbour = 0.0, top_neighbour = 0.0, this_node = 0.0;
- double left_coeff = 1.0, right_coeff = 1.0, bottom_coeff = 1.0, top_coeff = 1.0, this_coeff = 1.0;
- switch (IsNodeInternalCornerOrSide(glob_x, glob_y)){
- case 2:
- //left bottom corner
- neighbours["left"] = false;
- neighbours["bottom"] = false;
- break;
- case 4:
- //left top corner
- neighbours["left"] = false;
- neighbours["top"] = false;
- break;
- case 6:
- //right bottom corner
- neighbours["right"] = false;
- neighbours["bottom"] = false;
- break;
- case 8:
- //right top corner
- neighbours["right"] = false;
- neighbours["top"] = false;
- break;
- case 1:
- //left side
- neighbours["left"] = false;
- break;
- case 3:
- //right side
- neighbours["right"] = false;
- break;
- case 5:
- //bottom side
- neighbours["bottom"] = false;
- break;
- case 7:
- //top side
- neighbours["top"] = false;
- break;
- case 0:
- //internal
- break;
- default:
- cout << "[ERROR]: Bad global coords compute matrix. Global:" << glob_x << " " << glob_y<<endl;
- }
- if (!neighbours["left"]){
- right_coeff = 2.0;
- left_coeff = 0.0;
- }
- if (!neighbours["right"]){
- left_coeff = 2.0;
- right_coeff = 0.0;
- }
- if (!neighbours["bottom"]){
- top_coeff = 2.0;
- bottom_coeff = 0.0;
- }
- if (!neighbours["top"]){
- bottom_coeff = 2.0;
- top_coeff = 0.0;
- }
- if (neighbours["left"]){
- left_coeff *= -k(x_left + (glob_x - 0.5) * h_x, y_bottom + glob_y * h_y) / pow(h_x, 2);
- left_neighbour = Get(glob_x - 1, glob_y);
- }
- if (neighbours["right"]){
- right_coeff *= -k(x_left + (glob_x + 0.5) * h_x, y_bottom + glob_y * h_y) / pow(h_x, 2);
- right_neighbour = Get(glob_x + 1, glob_y);
- }
- if (neighbours["bottom"]){
- bottom_coeff *= -k(x_left + glob_x * h_x, y_bottom + (glob_y - 0.5) * h_y) / pow(h_y, 2);
- bottom_neighbour = Get(glob_x, glob_y - 1);
- }
- if (neighbours["top"]){
- top_coeff *= -k(x_left + glob_x * h_x, y_bottom + (glob_y + 0.5) * h_y) / pow(h_y, 2);
- top_neighbour = Get(glob_x, glob_y + 1);
- }
- this_coeff = q(x_left + glob_x * h_x, y_bottom + glob_y * h_y) - left_coeff - right_coeff - bottom_coeff - top_coeff;
- this_node = Get(glob_x, glob_y);
- result = left_coeff * left_neighbour +
- right_coeff * right_neighbour +
- bottom_coeff * bottom_neighbour +
- top_coeff * top_neighbour +
- this_coeff * this_node;
- return result;
- }
- double GetNodeFromB(int current_node_global_offset_x, int current_node_global_offset_y) {
- int glob_x = current_node_global_offset_x;
- int glob_y = current_node_global_offset_y;
- double result = 0.0;
- switch (IsNodeInternalCornerOrSide(glob_x, glob_y)){
- case 2:
- //left bottom corner
- result = F(x_left, y_bottom) + 2.0 / h_x * psi_L(y_bottom) + 2.0 / h_y * psi_B(x_left);
- break;
- case 4:
- //left top corner
- result = F(x_left, y_top) + 2.0 / h_x * psi_L(y_top) + 2.0 / h_y * psi_T(x_left);
- break;
- case 6:
- //right bottom corner
- result = F(x_right, y_bottom) + 2.0 / h_x * psi_R(y_bottom) + 2.0 / h_y * psi_B(x_right);
- break;
- case 8:
- //right top corner
- result = F(x_right, y_top) + 2.0 / h_x * psi_R(y_top) + 2.0 / h_y * psi_T(x_right);
- break;
- case 1:
- //left side
- result = F(x_left, y_bottom + glob_y * h_y) + 2.0 / h_x * psi_L(y_bottom + glob_y * h_y);
- break;
- case 3:
- //right side
- result = F(x_right, y_bottom + glob_y * h_y) + 2.0 / h_x * psi_R(y_bottom + glob_y * h_y);
- break;
- case 5:
- //bottom side
- result = F(x_left + glob_x * h_x, y_bottom) + 2.0 / h_y * psi_B(x_left + glob_x * h_x);
- break;
- case 7:
- //top side
- result = F(x_left + glob_x * h_x, y_top) + 2.0 / h_y * psi_T(x_left + glob_x * h_x);
- break;
- case 0:
- //internal
- result = F(x_left + glob_x * h_x, y_bottom + glob_y * h_y);
- break;
- default:
- cout << "[ERROR]: Bad global coords compute matrix. Global:" << glob_x << " " << glob_y <<endl;
- }
- return result;
- }
- double func_inner_product(vector<double> &v1, vector<double> &v2) {
- double prod = 0;
- // rho = 1 in case of interior points
- for (int i = 1; i < M; ++i)
- for (int j = 1; j < N; ++j)
- prod += v1[j * (M + 1) + i] * v2[j * (M + 1) + i];
- // rho = 0.5 in case of left and right sides
- for (int j = 1; j < N; ++j)
- prod += 0.5 * (v1[j * (M + 1)] * v2[j * (M + 1)] + v1[j * (M + 1) + M] * v2[j * (M + 1) + M]);
- // rho = 0.5 in case of top and bottom sides
- for (int i = 1; i < M; ++i)
- prod += 0.5 * (v1[N * (M + 1) + i] * v2[N * (M + 1) + i] + v1[i] * v2[i]);
- // rho = 0.25 in case of corner points
- prod += 0.25 * (v1[0] * v2[0] + v1[M] * v2[M] + v1[N * (M + 1)] * v2[N * (M + 1)]
- + v1[N * (M + 1) + M] * v2[N * (M + 1) + M]);
- prod *= h_x * h_y;
- return prod;
- }
- double norm(vector<double> &v) {
- return sqrt(func_inner_product(v, v));
- }
- // solving system by method of least residuals
- void DoIterations() {
- double tau = 0.0, eps = 1e-06;
- vector<double> r((M + 1) * (N + 1), 0.0), Ar((M + 1) * (N + 1), 0.0);
- vector<double> old_w;
- int iterations = 0;
- do {
- // filling residual
- //for (int i = 0; i < sz; ++i)
- for(int i = 0; i < M + 1; ++i){
- for(int j = 0; j < N + 1; ++j){
- //r[i] = inner_product(A[i], w) - B[i];
- r[j * (M + 1) + i] = ComputeMagicInnerProductA_iw(i, j) - GetNodeFromB(i, j);
- }
- }
- old_w = w;
- w = r;
- // filling Ar which is product of matrix A and vector r
- //for (int i = 0; i < sz; ++i)
- //Ar[i] = inner_product(A[i], r);
- for(int i = 0; i < M + 1; ++i){
- for(int j = 0; j < N + 1; ++j){
- Ar[j * (M + 1) + i] = ComputeMagicInnerProductA_iw(i, j);
- }
- }
- w = old_w;
- // computing iteration parameter
- tau = func_inner_product(Ar, r) / func_inner_product(Ar, Ar);
- // computing new iteration of w
- for (int i = 0; i < (M + 1) * (N + 1); ++i)
- w[i] -= tau * r[i];
- ++iterations;
- cout << fabs(tau) * norm(r) - eps << endl;
- } while (fabs(tau) * norm(r) >= eps);
- cout << "Number of iterations: " << iterations << endl;
- }
- // computing norm of difference between approximate and exact solutions
- double Compare(vector<double> &exact) {
- int sz = exact.size();
- // filling difference vector
- vector<double> diff(sz, 0.0);
- for (int i = 0; i < sz; ++i)
- diff[i] = w[i] - exact[i];
- return norm(diff);
- }
- int main(int argc, char* argv[]) {
- x_left = -2, x_right = 3;
- y_bottom = -1, y_top = 4;
- N = atoi(argv[1]);
- M = N;
- h_x = (x_right - x_left) * 1.0 / M;
- h_y = (y_top - y_bottom) * 1.0 / N;
- vector<double> exact((M + 1) * (N + 1), 0.0);
- CalculateExactSolution(exact);
- w.resize((M + 1) * (N + 1));
- DoIterations();
- cout << "[INFO]: Diff with exact solution:" << Compare(exact) << endl;
- return 0;
- }
Add Comment
Please, Sign In to add comment