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.

  • $\begingroup$ While your idea is very good and totally valid in mathematical sense, there are two problems: First, how close is the function to 0 does not really tell how close is the argument to the root (take 0.00001x = 0.00001, the error r you proposed is about 10^-6 even when x is within 10^-3 of 1 - the true solution). For this we would somehow need to divide by "the derivative of the linear operator", which I can't grasp at this moment, and certainly don't know how to implement it. $\endgroup$ – user16320 Sep 6 '19 at 21:05
  • $\begingroup$ Second problem: to compute r (or max or whatever), after each iteration I would need to sweep the grid again (since in the first sweep, the values stored in f are mixed: some are from the previous iteration, some are old, so for the current node, the equation Ax = b is fulfilled (with surrounding nodes in the current state of the iteration), but it's not fulfilled precisely, when the iteration is finished). The point is: another sweep of the grid essentially renders the iterations twice as slow. Computation of the error should not take a significant portion of the total time taken. $\endgroup$ – user16320 Sep 6 '19 at 21:05
  • $\begingroup$ I think there's a point of confusion here, as what I suggested is the proper way, and when d is equal to 0, the x vector is the solution of the linear system. Furthermore, the derivative of a linear operator with respect to the input is simply the linear operator applied to the derivative of the input-- because the operator is linear. $\endgroup$ – EMP Sep 6 '19 at 21:09
  • $\begingroup$ Okay, let's revisit my little example. My linear system is 1x1 and looks like: 0.00001x = 0.00001. For x = 1.01 (one percent within the correct solution), your method of computing d, that is 0.00001 - 0.00001x gives 10^-7. Ergo, the "error estimate" (10^-7) does not even closely correspond to the actual error (one percent). $\endgroup$ – user16320 Sep 6 '19 at 21:12
  • $\begingroup$ You don't compute the error after every step, you compute every hundred, or at a number equal to the number of DOFs of the system. You're using Gauss-Seidel, you don't need the error at every step because it provides very little information. Your current implementation isn't in fact an error estimate. $\endgroup$ – EMP Sep 6 '19 at 21:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.