How to set the boundary conditions of 3D Poisson Equation

조회 수: 30 (최근 30일)
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;
  1. Dirichlet boundary condition are applied at the top and bottom of the planes of the rectangular grid.
  2. Electric potential is to be incorporated by setting and , where h is the height of the simulation box.
  3. 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
David Martín Oliva Zavala
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:)

댓글을 달려면 로그인하십시오.

채택된 답변

Dinesh Yadav
Dinesh Yadav 2020년 3월 5일
편집: Dinesh Yadav 2020년 3월 5일
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
Dinesh Yadav
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
Anthony Owusu-Mensah 2020년 3월 5일
편집: Anthony Owusu-Mensah 2020년 3월 5일
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

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Graphics Object Programming에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by