Neumann boundary with finite differences

[JohanL]

I tried to solve laplace equation for the steady state temperature over a rectangle with both neumann and dirichlet boundary conditions.
For the part of the rectangle with neumann boundary condition(normal derivative = 0) i used

U(k,p)=(2*U(k+1,p)+U(k,p+1)+U(k,p-1))/4

Is this correct?
The result i get is different from when i use Matlab PDE toolbox...and i think I am doing something wrong with the neumann boundary condition.

 

[J77]

Assuming you used discretisation:

\frac{1}{h^2}(u_{i+1}-2u_i+u_{i-1})

Say your b.c. is at X=0 (T for temp. etc.):

\frac{\partial T(0,t)}{\partial X}=0

For the discretisation, you don't know u_{i-1} but you can rewrite the discretisation as:

\frac{1}{h}\left(\frac{u_{i+1}-u_i}{h}-\frac{u_i-u_{i-1}}{h}\right)

Where the \frac{u_i-u_{i-1}}{h} bit is given by the b.c.

 

[JohanL]

I did it like this

if given homogenous neumann boundary condition for the points (0,j)

Laplace equation gives

4*T(0,j)=T(1,j)+T(-1,j)+T(0,j+1)+T(0,j-1) (*)

then

dT/dx=0 for those boundary points and i used

dT/dx=[T(1,j)-T(-1,j)]/2*h

and (-1,j) is outside the domain so i use

T(-1,j)=T(1,j)-2*h*dT/dx

and substituting this into (*) gives

4*T(0,j)=T(1,j)+T(1,j)-2*h*dT/dx+T(0,j+1)+T(0,j-1)

but dT/dx=0 which gives

4*T(0,j)=2*T(1,j)+T(0,j+1)+T(0,j-1)

______________________________________
in the code i have T=U and (i,j)=(k,p). The neumann boundary is for k=1 or Nx and 12<p<20 and otherwise dirichlet boundary conditions.

Nmax=500; eps=10^(-3); n=0; delta=1; % parametrar
while abs(delta)>eps & n<Nmax
n=n+1
for k=2:Nx-1
for p=2:Ny-1

if (k==1) & (12<p<20)
dU(k,p)=(2*U(k+1,p)+U(k,p+1)+U(k,p-1))/4 ...
-U(k,p);
U(k,p)=U(k,p)+om*dU(k,p);

elseif k==Nx & 12<p<20
dU(k,p)=(2*U(k-1,p)+U(k,p+1)+U(k,p-1))/4 ...
-U(k,p);
U(k,p)=U(k,p)+om*dU(k,p);

elseif (k~=1) | (k~=Nx)

dU(k,p)=(U(k+1,p)+U(k-1,p)+U(k,p+1)+U(k,p-1))/4 ...
-U(k,p);
U(k,p)=U(k,p)+om*dU(k,p);

else
U(k,p)=U(k,p)
end
end
end
delta=max(max(dU)); disp(delta);
end

 

[Clausius2]

I would rather leave interior domain points to Laplace discretization, and in the Neumman boundary I would chose either a 2nd order central derivative approx centered just on the boundary point and applying an image point at -1, or a one-sided 2nd order derivative in which you obtain the boundary point value as a function only of the inner points.

 

[JohanL]

I would rather leave interior domain points to Laplace discretization, and in the Neumman boundary I would chose either a 2nd order central derivative approx centered just on the boundary point and applying an image point at -1.

Isnt that what I am doing? (if you by image point mean phantom point?)

Or is it something wrong with my approach?

In the code its of course is
for
k=1:Nx

First i calculate for the points on the boundary which have a neumann boundary condition:

if (k==1) & (12<p<20)
dU(k,p)=(2*U(k+1,p)+U(k,p+1)+U(k,p-1))/4 ...
-U(k,p);
U(k,p)=U(k,p)+om*dU(k,p);

elseif k==Nx & (12<p<20)
dU(k,p)=(2*U(k-1,p)+U(k,p+1)+U(k,p-1))/4 ...
-U(k,p);
U(k,p)=U(k,p)+om*dU(k,p);

where om is a relaxation factor to speed up the calculations.

Then i calculate for the interior points

elseif (k~=1) | (k~=Nx)

dU(k,p)=(U(k+1,p)+U(k-1,p)+U(k,p+1)+U(k,p-1))/4 ...
-U(k,p);
U(k,p)=U(k,p)+om*dU(k,p);

and finally

else
U(k,p)=U(k,p)

for those boundary points with dirichlet boundary conditions
Nx = 36; Ny = 32;

U=zeros(Nx,Ny);
U(1:Nx,Ny)=100;
U(1:Nx,1)=20;
U(1,20:32)=100;
U(Nx,20:32)=100
U(1,1:12)=20;
U(Nx,1:12)=20;