Solving Laplace equation using Finite difference

Greetings all,

I'm trying to solve the following problem using a finite differnce iterative scheme. I wrote the following code which seems to give me a solution that does not vary with changing the length. I end up getting the same solution regardless of the domain length which contradicts the solution from the pde tool. My code is below, your help is appericiated. I think the problem is in how I implement the loop maybe, I'm not sure.

clear all;

clc;

% Apply user input here

L = 3;

H = 1;

pts = 80;

% In case of undivisible by 3 grid

int_part = fix(pts/3);

rem_part = pts - (int_part*3);

% Grid and tolerance

tol = 1.0e-5;

if rem_part ~= 0

x1 = int_part;

x2 = (x1 * 2);

x3 = (x1 * 3) + rem_part;

else

x1 = int_part;

x2 = x1*2;

x3 = x1*3;

end

x = 0:L/(pts-1):L;

y = 0:H/(pts-1):H;

%[X,Y] = meshgrid(x,y);

%% Initialization

U = zeros(pts);

U0 = zeros(pts);

eta = 1; % dummy

%% boundary conditions

r1 = 2:x1;

r2 = (x1+1):x2;

r3 = (x2+1):(x3-1);

rx = 2:pts-1;

U(pts,1:x1) = -100;

U(1,r2) = 100;

U(pts,(x2+1):(x3)) = -100;

%% Solution loop

i=2:pts-1;

j=2:pts-1;

while eta > tol

U(i,j) = 0.25*(U(i-1,j)+U(i+1,j)+U(i,j-1)+U(i,j+1)); % Inner Domain

U(rx,1) = 0.25*(U(rx,1)+U(rx,1)+U(rx,2)+U(rx,2)); % b.c -> ux = 0

U(rx,pts) = 0.25*(U(rx,pts)+U(rx,pts)+U(rx,pts-1)+U(rx,pts-1)); % b.c -> ux = 0

U(1,r1) = 0.25*(U(2,r1)+U(2,r1)+U(1,r1-1)+U(1,r1+1)); % b.c -> uy = 0

U(1,r3) = 0.25*(U(2,r3)+U(2,r3)+U(1,r3-1)+U(1,r3+1)); % b.c -> uy = 0

U(pts,r2) = 0.25*(U(pts-1,r2)+U(pts-1,r2)+U(pts,r2-1)+U(pts,r2+1)); % b.c -> uy = 0

U(1,1) = 0.25*(U(2,1)+U(2,1)+U(1,2)+U(1,2)); % corner point

U(1,pts) = 0.25*(U(2,pts)+U(2,pts)+U(1,pts-1)+U(1,pts-1)); % corner point

eta = max(max(abs(U-U0)));

U0 = U;

end

contourf(x,y,U)

colorbar

colormap('jet')

Answer 1

[Torsten]

If L = 3 and H = 1, you cannot choose the same number of points in x and y direction if you update U as if you use a grid with equal cell sizes in x and y direction.

So either you choose nx = 3*ny for nx, ny being the number of cells (not points) in x and y direction, respectively, or you update your U values according to

(U(i+1,j)-2*U(i,j)+U(i-1,j))/dx^2 + (U(i,j+1)-2*U(i,j)+U(i,j-1))/dy^2 = 0

thus

U(i,j) = 1/(2/dx^2+2/dy^2) * ( (U(i+1,j)+U(i-1,j))/dx^2 + (U(i,j+1)+U(i,j-1))/dy^2 )

with

dx = L/nx, dy = H/ny

(Similar for the points where derivatives are prescribed).

And of course all the U's on the right-hand side of your updates

U(i,j) = 0.25*(U(i-1,j)+U(i+1,j)+U(i,j-1)+U(i,j+1));              % Inner Domain
U(rx,1) = 0.25*(U(rx,1)+U(rx,1)+U(rx,2)+U(rx,2));                 % b.c -> ux = 0
U(rx,pts) = 0.25*(U(rx,pts)+U(rx,pts)+U(rx,pts-1)+U(rx,pts-1));   % b.c -> ux = 0
U(1,r1) = 0.25*(U(2,r1)+U(2,r1)+U(1,r1-1)+U(1,r1+1));                % b.c -> uy = 0
U(1,r3) = 0.25*(U(2,r3)+U(2,r3)+U(1,r3-1)+U(1,r3+1));                % b.c -> uy = 0
U(pts,r2) = 0.25*(U(pts-1,r2)+U(pts-1,r2)+U(pts,r2-1)+U(pts,r2+1));  % b.c -> uy = 0  
U(1,1) = 0.25*(U(2,1)+U(2,1)+U(1,2)+U(1,2));                % corner point
U(1,pts) = 0.25*(U(2,pts)+U(2,pts)+U(1,pts-1)+U(1,pts-1));  % corner point

must be changed to U0's.

[Khaled Mohamed]

Thank you. That was the problem. I was confused. Your help is really apperciated.