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

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++)
{

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$$

• How did you treat your boundary conditions? Please write out your scheme at the boundary. – Paul Oct 18 '16 at 15:24
• What order are you computing? And are you using centered difference? – Charles Oct 18 '16 at 16:19
• This question is really unclear without an MCVE. Could you post one, please? – Kirill Oct 18 '16 at 16:38
• Can you also say on what domain you are solving? It would also help if you showed a convergence plot. What order do you get? – Wolfgang Bangerth Oct 18 '16 at 18:09
• @Paul:Please see edit Boundary conditions. – Gaurav Saxena Oct 18 '16 at 21:03

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):

• Actually, we use the continuous L2 norm. I think that's what you probably meant and compute. – Wolfgang Bangerth Oct 18 '16 at 22:11
• @user3302471: thanks. I misunderstood convergence. But I am still confused about the stopping criterion of the program. Whatever mesh size I take, the MAE = maximum absolute error = max( |analytical - calculated| ) refuses to fall below 0.314146. Shouldn't the absolute error keep decreasing ? Since I know my MAE never falls below 0.314146, will keeping the stopping criterion as the relative $l_{2}$ norm be useful ? – Gaurav Saxena Oct 19 '16 at 10:25
• The stopping criterion of your simulation should be on the residual of your system. In your case, you end up solving a problem of the form : A*x = b, where A is a matrix, x is the solution vector and b is the second member. When your iterative solver reaches convergence, this system will be solved close to the tolerance (that is between two iterations, the vector x will barely change). This means your simulation has reached it's tolerance (Your individual simulation has converged). You must not use the analytical solution to verify if a single simulation has reached convergence. – BlaB Oct 19 '16 at 14:41
• Just to add to this, I've found some papers by P.J. Roache very useful. I think one title is "QUANTIFICATION OF UNCERTAINTY IN COMPUTATIONAL FLUID DYNAMICS", and this author has several publications on this topic, all of which that I've read I've found very informative and useful. – Charles Oct 23 '16 at 3:41

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.

• as a side-note: for a non-uniform mesh, a good idea is to estimate the (numerical) convergence order from the relationship: error vs (total) degrees of freedom (DoFs) in your discretization. – GoHokies Oct 19 '16 at 19:34