$O(h^2)$ convergence for Elliptic PDE

Question

I am trying to solve an elliptic PDE in 2-D:

$$-\nabla^{2} u = 20tanh(10x-5)(10-10tanh(10x-5)^2) = f$$

I know that the solution is $u = tanh(10x-5)$ but I am unable to get $O(h^2)$ solution with a cell-centered, finite-difference on a uniform mesh with Dirichlet boundary conditions.

Discretization : central difference $$-[\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}+u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h^{2}}] = f_{i,j}$$ $$u_{i,j} = \frac{u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1}+f_{i,j}h^{2}}{4.0} $$

Boundary conditions : X-direction is from top to bottom, Y-direction is from left to right. Since I know the analytical solution : $u = tanh(10x-5)$, I substitute the value of $u$ for various discretised points at the boundary. Then at each update of $u$, I re-adjust the ghost cell as : $$\frac{u_{0,j}+u_{1,j}}{2.0} = tanh(10*0-5)$$ at the lower-X boundary (and similar expressions for high-X,low-Y, high-Y).

Plot: In 2-D/3-D I get a plot of the solution similar to the figure shown below. Please ignore the Adaptive Mesh.

Code: using Gauss-Seidel iterative scheme.

double update(double phi[][N+2],double rhs[][N+2], int Nx, int Ny, double dx, double kp )
{
int i,j;
double abs_err = 0.0, err, x; 

for(i = 1; i <= Nx-2; i++)
{
    x = (i-0.5)*dx ; 
    for(j = 1; j <= Ny-2; j++ )
    {
        phi[i][j] = (phi[i+1][j] + phi[i-1][j] + phi[i][j-1] + phi[i][j+1] + dx*dx*rhs[i][j])/4.0 ; 
        err = tanh(kp*(x-0.5)) - phi[i][j];
        if(fabs(err) > abs_err)
            abs_err = fabs(err); 
    }
}       
return abs_err ; 
}

void boundary(double phi[][N+2],int Nx,int Ny,double dx, double kp)
{
int i,j;
double x,y; 

//low-y boundary i.e. y = 0     
for(i = 1; i <= N+2-2; i++)
{
      x = (i-0.5)*dx ; 
      phi[i][0]      = 2 * tanh(kp*(x-0.5)) - phi[i][1]; 
}   
//high-y boundary i.e. y = 1    
for(i = 1; i <= N+2-2; i++)
{
    x = (i-0.5)*dx ; 
      phi[i][N+1]    = 2 * tanh(kp*(x-0.5)) - phi[i][N];        
}   
//low-x boundary i.e. x = 0
for(j = 1; j <= N+2-2; j++)
{        
      phi[0][j]    =     2 * tanh(kp*(0.0-0.5)) - phi[1][j];

}   
//high-x boundary i.e. x = 1    
for(j = 1; j <= N+2-2; j++)
{       
      phi[N+1][j]    =   2 * tanh(kp*(1.0-0.5)) - phi[N][j]; 

}   
}

void main()
{
double phi[N+2][N+2] ;                  //solution array with ghost boundaries
double rhs[N+2][N+2];                   //RHS i.e. -del^{2}u = RHS

double prob_x[2] = {0.0, 1.0};                  //X-direction problem domain
double prob_y[2] = {0.0, 1.0};                  //Y-axis problem domain
double x,y, kp = 10.0, abs_err; 

double dx, dy ;                                 //Mesh spacing
dx = dy = (prob_x[1] - prob_x[0])/N ; 

int i,j,ctr; 

for(i = 0; i <= N+2-1; i++)
    for(j=0; j <= N+2-1; j++)
        {
            phi[i][j] = 0.0;
            rhs[i][j] = 0.0; 
        }

for(i = 1; i <= N; i++)
{
    x = (i-0.5)*dx ; 
    for(j = 1; j<= N; j++)
    rhs[i][j] = 20 * tanh(10.0*x-5.0) * (10.0 - 10.0 * tanh( (10.0*x-5.0) * (10.0*x-5.0) ) ) ; 
}

boundary(phi,N+2,N+2,dx,kp);

ctr = 1; 
abs_err = 1.0; 

//while(abs_err > dx*dx)
for(ctr = 1; ctr <= 50000; ctr++)
{

    abs_err = update(phi,rhs,N+2,N+2,dx,kp);                //Updates phi 
    boundary(phi,N+2,N+2,dx,kp); 
    //ctr = ctr+1; 
}   
}

Convergence : I don't know the order of convergence I am getting. Sorry I don't know how to calculate that.

I do get solutions close to the analytic solutions but not the order $O(h^2)$ ($h$ = mesh points spacing) convergence. Can anyone shed some light on this ?

EDITS again Convergence Rate: Thanks for all the answers and helpful comments. I defined $r_{0} = ||\bar{x} - x^{(0)}||_{2}$ and $r_{k} = ||\bar{x} - x^{(k)}||_{2}$ where $\bar{x}$ is the actual analytical solution and $x^{k}$ is calculated solution at the $k^{th}$ iteration. I agree that properly the stopping criterion should be the relative norm of residuals but here I took the norm of error measures from the analytical solution. I am finally getting these values.

• Mesh size : 4x4, $\frac{r_{k}}{r_{0}} = 0.31$
• Mesh size : 8x8, $\frac{r_{k}}{r_{0}} = 0.092459$
• Mesh size : 16x16, $\frac{r_{k}}{r_{0}} = 0.01122066$
• Mesh size : 32x32, $\frac{r_{k}}{r_{0}} = 0.0026179$
• Mesh size : 64x64, $\frac{r_{k}}{r_{0}} = 0.000646645$
• Mesh size : 128x128, $\frac{r_{k}}{r_{0}} = 0.00016118577$

I think this qualifies for a quadratic convergence.

Mistake : I made a mistake in the expression for RHS which I had obtained from doubly differentiating $u = tanh(10x-5)$. I used Maple for it (first time user) and I did not know that Maple writes $tanh^{2}(10x-5)$ as $tanh(10x-5)^{2}$. Hence the correct PDE is :

$$-\nabla^{2} u = 20tanh(10x-5)(10-10tanh(10x-5)tanh(10x-5)) = f$$

Answer 1

You appear (and correct me if I am wrong, i do not wish to presume) to be confused on the notion of what is an order of convergence.

The order of convergence is the order at which your approximate solution (in this case, your simulation) approaches the theoretical solution.

In your case, the theoretical order of convergence of the centered scheme is 2. This means that if you refine your mesh by a factor of two, the error on the approximate solution, which is O(h^2) will decrease by a factor of 4.

Now, defining the error is not trivial. Generally we use the L2 norm for the error (http://mathworld.wolfram.com/L2-Norm.html) for continous problems (for a discontinous problem such as a shock wave, the L1 or Linf norm are also often used).

Concretly, to show that you get second order convergence, you need to calculate the L2 error between your numerical solution and the analytical solution for a number of meshes (say 10x10x10, 20x20x20, 30x30x30, etc.) and then carry out a linear regression such that (log(error) = slope * log (dx)) where dx is the mesh spacing. The slope will be the order of convergence of your implemented scheme, and in your case should be 2.

Such a graph should look like this (taken from a paper by me):

Answer 2

I think if you are interested to find the order of convergence then you should vary the mesh spacing ($h$). Now as you vary the mesh spacing you can obtain the error of the numerical solution. Then if you do a log plot of the error and $h$, you should get a straight line with slope $2$ I guess. And that is how you can estimate the order of convergence, quadratic in this case I presume.