#include #include #include class Node { private: double v; int fixed; public: Node() : v(0.0), fixed(0) { // empty } Node(double v, int fixed) : v(v), fixed(fixed) { // empty } double getV() const { return v; } void setV(double nv) { v = nv; } int getFixed() const { return fixed; } void setFixed(int nf) { fixed = nf; } }; void setBoundary(std::vector>& m) { m[1][1].setV(1.0); m[1][1].setFixed(1); m[6][7].setV(-1.0); m[6][7].setFixed(-1); } double calculateDifference(const std::vector>& m, std::vector>& d, const int w, const int h) { double total = 0.0; for (int i = 0; i < h; ++i) { for (int j = 0; j < w; ++j) { double v = 0.0; int n = 0; if (i > 0) { v += m[i - 1][j].getV(); n++; } if (j > 0) { v += m[i][j - 1].getV(); n++; } if (i + 1 < h) { v += m[i + 1][j].getV(); n++; } if (j + 1 < w) { v += m[i][j + 1].getV(); n++; } v = m[i][j].getV() - v / n; d[i][j].setV(v); if (m[i][j].getFixed() == 0) { total += v * v; } } } return total; } double iter(std::vector>& m, const int w, const int h) { using namespace std; vector> d; for (int i = 0; i < h; ++i) { vector t(w); d.push_back(t); } double curr[] = { 0.0, 0.0, 0.0 }; double diff = 1e10; while (diff > 1e-24) { setBoundary(m); diff = calculateDifference(m, d, w, h); for (int i = 0; i < h; ++i) { for (int j = 0; j < w; ++j) { m[i][j].setV(m[i][j].getV() - d[i][j].getV()); } } } for (int i = 0; i < h; ++i) { for (int j = 0; j < w; ++j) { int k = 0; if (i != 0) ++k; if (j != 0) ++k; if (i < h - 1) ++k; if (j < w - 1) ++k; curr[m[i][j].getFixed() + 1] += d[i][j].getV()*k; } } return (curr[2] - curr[0]) / 2.0; } const int S = 10; int main() { using namespace std; vector> mesh; for (int i = 0; i < S; ++i) { vector t(S); mesh.push_back(t); } double r = 2.0 / iter(mesh, S, S); cout << "R = " << setprecision(15) << r << '\n'; return 0; }