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