# Jacobi iteration for finite difference: when to stop?

I implemented a finite difference scheme to solve Poisson's equation in a 2D grid in C. I solve the system by using Jacobi iteration. Everything works fine until I use a while loop to check whether it is time to stop iterating or not (with for loops is easy). On the notes I am following there is written that I have to compute the following:

$$\delta = max ||uNew - u||,$$

for $1 <= i,j <= n$. With $uNew$ being the current solution and $u$ the previous iteration. Obviously they are 2D arrays.

I tried to do the following:

// Iterations
double delta=(tau+1), temp_delta;
#ifdef _OPENMP
wt1=omp_get_wtime();
#endif
do {

for(i=1;i<y-1;i++) {
for(j=1;j<x-1;j++) {
uNew[i*x+j] = 0.25 * (u[(i-1)*x+j] + u[i*x+(j+1)] + u[i*x+(j-1)] + u[(i+1)*x+j] - dx*dy*func(i,j,dx,dy));
}
}

// Boundary conditions using g(x,y)
for(j=0;j<x;j++) {
uNew[j] = gunc(0,j,dx,dy);
uNew[(y-1)*x+j] = gunc(y-1,j,dx,dy);
}
for(i=0;i<y;i++) {
uNew[i*x] = gunc(i,0,dx,dy);
uNew[i*x+(x-1)] = gunc(i,x-1,dx,dy);
}

// Check if to terminate Jacobi iteration
for(i=1;i<y-1;i++) {
for(j=1;j<x-1;j++) {
temp_delta = abs(uNew[i*x+j]-u[i*x+j]);
printf("%f ", temp_delta);
if (delta <= temp_delta) {
delta = temp_delta;
}
}
}

// Update solution
for(i=0;i<y;i++) {
for(j=0;j<x;j++) {
u[i*x+j] = uNew[i*x+j];
}
}
} while(delta > tau);


where $\tau$ is the tolerance, $\delta$ is the result of the above formula, $temp\_delta$ is used to find the maximum and $uNew$ and $u$ are just the matrices containing the solutions at the grid points. The problem is that it's not working. Can somebody give me a hint, please?

• What does the output of your printf statement say, compared to tau? And what is not working? Convergence or something else? – Bort Nov 21 '16 at 9:35
• @Bort It stops but at the wrong iteration since the solution is not correct. – wrong_path Nov 21 '16 at 9:41

For solving systems you shouldn't be comparing the results between each iteration but rather computing the residual.

If you consider the matrix representation to be in the standard form:

$A\cdot y = b$

Then you can define the residual at some iteration ($y_{i}$)

$r =b - A\cdot y_i$

Then you 'just' need to determine a stopping method based on the residual vector, some choices are

$\delta < max(|r|)$

$\delta < sum(|r|)$

Wikipedia article on the residual

• Thanks for the explanation but I do not have a system in this case. – wrong_path Nov 21 '16 at 12:14
• In fact, you have a system of linear algebraic equations, but you are solving it in a matrix free form. Your residual is a vector where each component is equal u[ix+j] - 0.25 * (u[(i-1)*x+j] + u[ix+(j+1)] + u[ix+(j-1)] + u[(i+1)*x+j] - dxdy*func(i,j,dx,dy)); – Peter Frolkovič Nov 21 '16 at 13:54
• Another common criteria is to ensure $r_i/r_0 \lt tolerance$, where $r_i$, $r_0$ are the ith and initial residuals. This is sometimes used to estimate solutions at each time step of some time marching method. – Charles Dec 6 '16 at 6:16