równanie Poissona

0

Czy byłby ktoś tak uprzejmy i wytłumaczył mi po kolei co ten program robi?

// Poisson's equation in 2-D using Jacobi, Gauss-Seidel,
//     or Successive Over Relaxation

#include <cmath>
#include <cstdlib>
#include <ctime>
#include <iostream>
#include <fstream>
#include <vector>
using namespace std;

int L;                        // number of interior points in x and y
vector< vector<double> > psi; // potential to be found
vector< vector<double> > rho; // given charge density

double h;                     // lattice spacing
vector< vector<double> > psiNew;   // new potential after each step
int steps;                    // number of iteration steps
double accuracy;              // desired accuracy in solution
double omega;                 // overrelaxation parameter

void initialize() {

    // create matrices with zero entries
    int N = L + 2;            // interior points + 2 boundary points
    psi = psiNew = rho
      = vector< vector<double> >(N, vector<double>(N, 0.0));

    h = 1 / double(L + 1);    // assume physical size in x and y = 1
    double q = 10;            // point charge
    int i = L / 2;            // center of lattice
    rho[i][i] = q / (h * h);  // charge density

    steps = 0;
}

void Jacobi() {

    // Jacobi algorithm for a single iterative step
    for (int i = 1; i <= L; i++)
    for (int j = 1; j <= L; j++)
        psiNew[i][j] = 0.25 * (psi[i - 1][j] + psi[i + 1][j] +
                               psi[i][j - 1] + psi[i][j + 1] +
                               h * h * rho[i][j]);
}

double relativeError() {

    double error = 0;         // average relative error per lattice point
    int n = 0;                // number of non-zero differences

    for (int i = 1; i <= L; i++)
    for (int j = 1; j <= L; j++) {
        if (psiNew[i][j] != 0)
        if (psiNew[i][j] != psi[i][j]) {
            error += abs(1 - psi[i][j] / psiNew[i][j]);
            ++n;
        }
    }
    if (n != 0)
        error /= n;

    return error;
}

void GaussSeidel() {

    // copy psi to psiNew
    for (int i = 1; i <= L; i++)
    for (int j = 1; j <= L; j++)
        psiNew[i][j] = psi[i][j];

    // Gauss-Seidel update
    for (int i = 1; i <= L; i++)
    for (int j = 1; j <= L; j++)
        psiNew[i][j] = 0.25 * (psiNew[i - 1][j] + psiNew[i + 1][j] +
                               psiNew[i][j - 1] + psiNew[i][j + 1] +
                               h * h * rho[i][j]);
}

void SuccessiveOverRelaxation() {

    // update even sites
    for (int i = 1; i <= L; i++)
    for (int j = 1; j <= L; j++)
    if ((i + j) % 2 == 0)
        psiNew[i][j] = (1 - omega) * psi[i][j] + omega / 4 *
                       (psi[i - 1][j] + psi[i + 1][j] +
                        psi[i][j - 1] + psi[i][j + 1] +
                        h * h * rho[i][j]);

    // update odd sites
    for (int i = 1; i <= L; i++)
    for (int j = 1; j <= L; j++)
    if ((i + j) % 2 != 0)
        psiNew[i][j] = (1 - omega) * psi[i][j] + omega / 4 *
                       (psiNew[i - 1][j] + psiNew[i + 1][j] +
                        psiNew[i][j - 1] + psiNew[i][j + 1] +
                        h * h * rho[i][j]);
}

void iterate(void (*method)()) {

    while (true) {
        method();
        ++steps;
        double error = relativeError();
        if (error < accuracy)
            break;
        vector< vector<double> > swap = psi;
        psi = psiNew;
        psiNew = swap;
    }
    cout << " Number of steps = " << steps << endl;
}

int main() {

    cout << " Iterative solution of Poisson's equation\n"
         << " ----------------------------------------\n";
    cout << " Enter number of interior points in x or y: ";
    cin >> L;

    initialize();

    cout << " Enter desired accuracy in solution: ";
    cin >> accuracy;

    clock_t t0 = clock();
    cout << " Enter 1 for Jacobi, 2 for Gauss Seidel, 3 for SOR: ";
    int choice;
    cin >> choice;
    switch (choice) {
    case 1:
        iterate(Jacobi);
        break;
    case 2:
        iterate(GaussSeidel);
        break;
    case 3:
        cout << " Enter overrelaxation parameter omega: ";
        cin >> omega;
        iterate(SuccessiveOverRelaxation);
        break;
    default:
        cout << " Jacobi: " << endl;
        iterate(Jacobi);
        cout << " Gauss-Seidel: " << endl;
        initialize();
        iterate(GaussSeidel);
        omega = 2 / (1 + 4 * atan(1.0) / double(L));
        cout << " Successive Over Relaxation with theoretical optimum omega = "
             << omega << endl;
        initialize();
        iterate(SuccessiveOverRelaxation);
        break;
    }

    clock_t t1 = clock();
    cout << " CPU time = " << double(t1 - t0) / CLOCKS_PER_SEC
         << " sec" << endl;

    // write potential to file
    cout << " Potential in file poisson.data" << endl;
    ofstream file("poisson.data");
    for (int i = 0; i < L + 2; i++) {
        double x = i * h;
        for (int j = 0; j < L + 2; j++) {
            double y = j * h;
            file << x << '\t' << y << '\t' << psi[i][j] << '\n';
        }
        file << '\n';
    }
    file.close();
}
 
0

moja (i sądzę ze większość forumowiczów) taka uprzejmość kosztuję. Dział praca z podaną kwotą zaprasza

0

kod nie jest napisany chaotycznie (tylkio main jest przydługawe) i zawiera odpowiednią ilość komentarzy, nić tylko googlać co ważniejsze komentarze i wszystko powinno być jasne.

1 użytkowników online, w tym zalogowanych: 0, gości: 1