# The proper way to assess the error of Jacobi iteration (for 2D Poisson equation)?

Motivation: I'm using 2D regular grid (it's actually a quadtree but I can still treat it as a finite difference thing if I weight-average the solution over smaller scale cells for the purpose of estimating Laplacian in the neighbor points) to discretize my Poisson equation and I'd like to solve it iteratively. The discretize Poisson equation $$\Delta f = g$$ reads

$$\frac{f_{i+1,j} + f_{i,j+1} + f_{i,j+1} + f_{i,j-1} - 4f_{i,j}}{h^2} = g_{i,j}$$ or, as an iterative method $$f_{i,j} = \frac{1}{4} \left( h^2 g_{i,j} - f_{i+1,j} - f_{i,j+1} - f_{i,j+1} - f_{i,j-1} \right)$$ meaning, that if $$g_{i,j}$$ are given, I traverse the grid and update $$f_{i,j}$$ based on the neighbour values of $$f$$ from the previous iteration (and since I don't keep the old values in memory, and always use the most recent values, I think the method is actually called "Gauss-Siedel")

I believe that this iterative method is identical to the Jacobi/Gauss-Siedel method used to solve linear equations (when formulated via matrix and RHS vector).

Error estimation: There are several formulas to assess the error (let's call it $$d$$) of the solution. All have the same property:

solution converges $$\implies$$ $$d$$ goes to zero

but each has some quirks that I will now describe:

The first few formulas probably anyone tries is ($$N$$ is number of grid points, occasionally appearing $$\varepsilon$$ is to make sure denominators are not zero and $$\alpha$$ is some chosen exponent):

$$d = \sum_{i,j} \left| f_{i,j} - f^{(old)}_{i,j} \right|^\alpha$$ $$d = \sum_{i,j} \left| \frac{f_{i,j} - f^{(old)}_{i,j}}{|f_{i,j}| + \varepsilon} \right|^\alpha$$ $$d = \frac{1}{N} \sum_{i,j} \left| f_{i,j} - f^{(old)}_{i,j} \right|^\alpha$$

However, the problem of this is that it smooths out some possible local non-convergence (if one point or group of points is not converging as quickly as others, it will get killed with the $$1/N$$ term or something similar).

Hence, I started to use something like the following $$d = \max_{\substack{i, j}} \{ \left| f_{i,j} - f^{(\text{old})}_{i,j} \right| \}$$ (again, with possible variations like terms divided by the value of $$f_{i,j}$$ to make it dimensionless etc.)

The problem: even when the cutoff on $$d$$ is small, like, 0.01 (so the computation stops when $$d < 0.01$$), the solution still may be far from the true solution (certainly farther than within one percent), so it seems like the difference between successive iterations is not enough to truthfully assess the error between the iterative solution and the true solution. Every paper I read on this do not address the question "when to stop the iteration", it is somehow generally understood that the answer is "when it stops changing", but that might not be enough (sometimes if I let it run for twice as many steps, I still get a much better solution, but $$d$$ is already ridiculously small, like $$10^{-6}$$). Of course, I know it's hard to estimate how far we are from the solution, especially since knowing the true values of the solution would eliminate the need to solve the Poisson equation in the first place.

The question: What is the best upper bound $$d$$ on the error of Jacobi iterations (given we don't know the true solution, only these iterative solutions) that still satisfies:

iterative solution converges to the true solution $$\iff$$ $$d \to 0$$?

(and reflects the nature of the "closeness" of the approximate solution to the true solution more faithfully? For example, if $$f_{i,j}$$ is 5 percent off from the true solution, the value of $$d$$ will be something more realistic, like $$0.05$$, rather than $$0.001$$)

Based off what you said, you have a linear poisson equation, and you want a measure of how well your answer that you compute from the linear system corresponds to the solution of the finite difference equation, which is itself a linear system. Since you have a linear system that you are solving through a Gauss-Seidel iteration, you can look at the problem you're solving as $$Ax = b$$ and you can define a residual $$r = b - Ax$$ or with your notation: $$d = g - Af$$ When you've satisfied the system this d will be equal to 0, just like you want. This is easy to compute and gives you the actual error as opposed to your current suggestions, which while they are nice will have problems dealing with stiff systems that don't converge well.