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();
}