How to set the boundary conditions of 3D Poisson Equation

Question

Anthony Owusu-Mensah 2020-3-2

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/508604-how-to-set-the-boundary-conditions-of-3d-poisson-equation

评论： David Martín Oliva Zavala 2020-11-12

I am trying to compute the electric potential at point (x,y,z) by solving the 3D Poisson equation below using finite difference method.

I have the charge densities at various positions of (x,y,z).

Below are the boundary condtions;

Dirichlet boundary condition are applied at the top and bottom of the planes of the rectangular grid.
Electric potential is to be incorporated by setting and , where h is the height of the simulation box.
Neumann boundary conditions are also enforced at the remaining box interfaces by setting at faces with constant x , at faces with constant y, and at faces with constant z.

I did the implementation of the boundary conditions with this code and I would want to ascertain if

the implementation of the boundary condtions is right.

x1 = linspace(0,10,20);
y1 = linspace(0,10,20);
z1 = linspace(0,10,20);
V = zeros(length(x1),length(y1),length(z1)); 
%Dirichlet Boundary Conditions
%Top plane
V(:,end,end) = 0;
V(end,:,end) = 0;
V(end,end,:) = 0;
% Bottom plane
V(:,1,1) = 0;
V(1,:,1) = 0;
V(1,1,:) = 0;
%Incoporated electric potential
V(:,:,1)  = 0; 
V(:,:,end) = -40*z1(end);
%Neumann Boundary Condition
i  = 2:length(x1)-1;
j = 2:length(y1)-1;
k = 2:length(z1)-1;
V(i+1,j,k) = V(i-1,j,k);
V(i,j+1,k) = V(i,j-1,k);
V(i,j,k+1) = V(i,j,k-1);          
        

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

David Martín Oliva Zavala 2020-11-12

Hey Anthony, have you finished solving poisson eqn in 3D? I´m just starting this project and I may need some help in the future jeje:)

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Dinesh Yadav 2020-3-5

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/508604-how-to-set-the-boundary-conditions-of-3d-poisson-equation#answer_418643

编辑：Dinesh Yadav 2020-3-5

在 MATLAB Online 中打开

Hi Anthony I think you have wrongly defined the top and bottom planes.A plane is a 2-D sheet structure however the line of code 
V(:,end,end);

represent a single row vector, similar to a line and not a plane. As V is a 20x20x20 3-D matrix, lets assume the first plane facing us is top plane i.e (assuming x direction increases no of columns and y direction increases no of rows, i.e origin is at (20,1,1))

V(:,:,1);    % top plane        face with constant z
V(:,:,end);   % bottom plane    face with constant z
V(:,1,:);    % face with constant x
V(:,end,:)   % face with constant x
V(1,:,:);    % face with constant y
V(end,:,:)   % face with constant y
dx  = 2:length(x1)-1;
dy = 2:length(y1)-1;
dz = 2:length(z1)-1;
V(dy,dx+1,dz) = V(dy,dx-1,dz);   % dv/dx = 0
V(dy+1,dx,dz) = V(dy-1,dx,dz);   % dv/dy = 0
V(dy,dx,dz+1) = V(dy,dx,dz-1);   % dv/dz = 0

3 个评论
显示 1更早的评论隐藏 1更早的评论

Dinesh Yadav 2020-3-5

Yes, replace i with dy j with dx and k with dz. Just to make it more readable. I have corrected the above code.

Anthony Owusu-Mensah 2020-3-5

编辑：Anthony Owusu-Mensah 2020-3-5

在 MATLAB Online 中打开

Thank you very much. But I guess changing the position of the x and y indices won't affect the solution because I intend on maitaining the indeces for x and y when solving for the internal nodes else per your submission I would have to also change the x and y indices for that too. Below is how I intend to slove for the internal nodes

l  = 2:length(x1)-1;
m = 2:length(y1)-1;
n = 2:length(z1)-1;
for t = 1:1000
    Vold = V;
    %%Solve for the internal nodes
    for i = 2:length(x1)-1 
        for j  = 2:length(y1)-1 
            for k = 2:length(z1)-1 
            V(i,j,k) = 1/6*(Vold(i+1,j,k)+ Vold(i-1,j,k)+...
               Vold(i,j+1,k)+ Vold(i,j-1,k)+...
               Vold(i,j,k+1)+ Vold(i,j,k-1));
            end
        end
    end
    %Dirichlet Boundary Condtions
    %Neumann Boundary Conditions
    V(l+1,m,n) = V(l-1,m,n);
    V(l,m+1,n) = V(l,m-1,n);
    V(l,m,n+1) = V(l,m,n-1); 
end