steady state 2D Poisson's equation - boundary conditions

Hello,

I am trying to solve steady state 2D Poisson's equation for a two spreading ridges problem. I have introduced boundary conditions tto. However though I should not have a flow across the top and bottom boundaries as you can see in the resultant image i do get flow across those boundaries. This gives me bulls eye contours on the spreading ridges but I should be getting horse shoe shapes instead.

This code solves for the 2D steady-state Poisson’s equation using finite difference method. The Poisson’s equation is given by:

∇²phi(x,y)=Gstar(x,y)

where phi(x,y) is the unknown function and Gstar(x,y) is a known function. The solution is obtained on a square domain with Dirichlet boundary conditions. The potential should not have this bull's eye shape. Potential about the ridges should be low and increase towards the boundaries. At the top and the bottom boundaries near ridges I know I should get a shape similar to horse shoe but as you can see something is not correct in my boundary conditions. For the top and bottom boundaries dphi/dx should be 0, which I thought I am doing but I cannot get this to work.

I am not sure if this is a math problem or Matlab problem anymore. Please help me with solving this problem.

%% Define grid parameters
x = 0:5:500;    dx=x(2)-x(1);   nx=length(x);
y = 0:5:500;    dy=y(2)-y(1);   ny=length(y);

% Define ridge locations
r1_x = 200;
r1_y = 250:500;
r2_x = 400;
r2_y = 0:250;

% Create meshgrid for plotting
[xx,yy] = meshgrid(x,y);

upper_half = yy >= 250;
lower_half = yy < 250;

% Apply age values to upper and lower halves of grid
t_upper = abs(xx - r1_x)/4;
t_lower = abs(xx - r2_x)/4;

t = zeros(size(xx));
t(upper_half) = t_upper(upper_half);
t(lower_half) = t_lower(lower_half);
t=max(t,0.1);

% Compute derivative of g
G = 0.11./(2*sqrt(t));

% Ignore infinity VALUES
count=0;
Gsum=0; 

for i=1:nx
    for j=1:ny
        if G(i,j)<=0.11
            Gsum=Gsum+G(i,j);
            count=count+1;
        end
    end
end

% Calculate the mean
Gmean=Gsum/count;
Gstar=G-Gmean;

%% Compute potentials using Poisson's equation (del^2*potential = Gstar)

phi = zeros(nx,ny); % Initialize solution matrix

% Top boundary condition
for j=2:nx-1
      phi(1,j) =phi(1,j-1)+phi(1,j+1)+phi(2,j)+phi(2,j)-(dx^2*Gstar(1,j)/4);
end

% Bottom boundary condition
for j=2:nx-1
    phi(ny,j) =phi(ny,j-1)+phi(ny,j+1)+phi(ny-1,j)+phi(ny-1,j)-(dx^2*Gstar(ny,j)/4);
end

% left boundary condition
for i=2:ny-1
      phi(i,1) =phi(i,2)+phi(i,2)+phi(i-1,1)+phi(i+1,1)-(dy^2*Gstar(i,1)/4);
end

% right boundary condition
for i=2:ny-1
      phi(i,end) =phi(i,end-1)+phi(i,end-1)+phi(i-1,end)+phi(i+1,end)-(dy^2*Gstar(i,end)/4);
end

% Top right corner

phi(1,1) = (dx^2*dy^2/(2*(dx^2 + dy^2))) * ((2*(phi(2,1))/dx^2)+ (2*(phi(1,2))/dy^2) - Gstar(1,1)) ;

% Top left corners

phi(1,nx) = (dx^2*dy^2/(2*(dx^2 + dy^2))) * ((2*(phi(2,nx))/dx^2)+ (2*(phi(1,nx-1))/dy^2) - Gstar(1,nx)) ;

% Bottom right corner

phi(ny,1) = (dx^2*dy^2/(2*(dx^2 + dy^2))) * ((2*(phi(ny-1,1))/dx^2)+ (2*(phi(ny,2))/dy^2) - Gstar(ny,1));

% Bottom left corner

phi(ny,nx) = (dx^2*dy^2/(2*(dx^2 + dy^2))) * ((2*(phi(ny-1,nx))/dx^2)+ (2*(phi(ny,nx-1))/dy^2) - Gstar(ny,nx));

%% 

% Jacobi iteration method
maxiter = 10000; % maximum number of iterations
tol = 1e-6; % tolerance for convergence
for iter = 1:maxiter
    phi_old = phi;
    for i = 2:ny-1
        for j = 2:nx-1
            phi(i,j) = (dx^2*dy^2/(2*(dx^2 + dy^2))) * (((phi(i+1,j) + phi(i-1,j))/dx^2)+ ((phi(i,j+1) + phi(i,j-1))/dy^2) - Gstar(i,j)) ;
        end
    end
    E = max(abs(phi(:)-phi_old(:)));
    if E < tol % convergence check
        break;
    end
end

%% Convert potential into velocities (u and v components)
p = reshape(phi,ny,nx);
u = 0*p;
v = 0*p;

u(:,2:end-1) = -(p(:,3:end) - p(:,1:end-2))/(2*dx);
v(2:end-1,:) = (p(3:end,:) - p(1:end-2,:))/(2*dy);

scale=5;

figure; hold on;axis fill;
q=quiver(x,y,u,v,scale);
contour(x,y,p,LineWidth=1)

[Torsten]

You want to encourage us to study your question ? We are lazy people who don't like to deduce the underlying problem from your code. So you should write down what you are trying to solve (equation, boundary conditions, method of solutions).

It seems that your boundary values depend on the values in the interior of the domain. Thus you have to adapt them during the Jacobi iteration. And shouldn't the right-hand side of the update

phi(i,j) = (dx^2*dy^2/(2*(dx^2 + dy^2))) * (((phi(i+1,j) + phi(i-1,j))/dx^2)+ ((phi(i,j+1) + phi(i,j-1))/dy^2) - Gstar(i,j)) ;

be written in phi_old and not in phi ?

[Anitha Limann]

Thank you for your kind answer. When i include boundary conditions in to jacobi interation method I get empty plots.

[Torsten]

When i include boundary conditions in to jacobi interation method I get empty plots.

Most probably your Jacobi iterations diverge and you get NaN values for p which then are not plotted. But getting no plot is not a reason not to code correctly.

And what about changing phi to phi_old for the right-hand side ?

And what about a problem description ?