1
$\begingroup$

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?

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

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

$\endgroup$
  • $\begingroup$ Thanks for the explanation but I do not have a system in this case. $\endgroup$ – wrong_path Nov 21 '16 at 12:14
  • 4
    $\begingroup$ 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)); $\endgroup$ – Peter Frolkovič Nov 21 '16 at 13:54
  • 2
    $\begingroup$ 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. $\endgroup$ – Charles Dec 6 '16 at 6:16

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.