Numerical Integration over a surface, for the verification of gauss law

8 次查看（过去 30 天）
Georgios Demetriou
回答： SAI SRUJAN 2024-1-31
I am trying to Numerically integrate over a surface to verify Gauss Law, using Gauss Seidel algorithm. My method consists of the following code.
function charge=compute_charge(xi,xf,yi,yf,X,Y,Ex,Ey)
% Prendre une surface rectangulaire passant par
% les points de maillage et entourant une électrode.
% Utiliser la formule des trapèzes pour les intégrales.
dx=(X(2)-X(1));
dy=(Y(2)-Y(1));
xmin=floor(xi./dx);
xmax=ceil(xf./dx)+1;
ymin=floor(yi./dy);
ymax=ceil(yf./dy)+1;
integral = 0.;
for i=xmin:xmax
integral=integral+dx*(Ey(i,ymax)-Ey(i,ymin));% it can be changed if there is no need for 1/2
end
%charges sur les surfaces verticales
for i=ymin:ymax
integral=integral+dy*(Ex(xmax,i)-Ex(xmin,i));
end
epsilon0 = 8.85418782e-12;
charge = integral * epsilon0 * 1e12;
% TODO: intégrer le flux du champ électrique
fprintf('Q = %0.3f pC\n\n',epsilon0*integral*1e12)
end
However I tend to beleive that this code has certain errors.
Could you help me out.

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

回答（1 个）

SAI SRUJAN
SAI SRUJAN 2024-1-31
Hi Georgios,
I understand that you are trying to numerically integrate over a surface using Gauss Seidel algorithm.
Please refer to the following code outline to proceed further,
x_min = 0;
x_max = 1;
y_min = 0;
y_max = 2;
% Number of grid points
Nx = 100;
Ny = 200;
dx = (x_max - x_min) / (Nx - 1);
dy = (y_max - y_min) / (Ny - 1);
phi = zeros(Ny, Nx);
% Set the boundary conditions
phi(:, 1) = 1;
phi(:, end) = 0;
phi(1, :) = 0;
phi(end, :) = 0;
% Perform Gauss-Seidel iteration
max_iter = 1000;
tolerance = 1e-6;
for iter = 1:max_iter
phi_old = phi;
for i = 2:Nx-1
for j = 2:Ny-1
phi(j, i) = 0.25 * (phi(j, i-1) + phi(j, i+1) + phi(j-1, i) + phi(j+1, i));
end
end
if max(abs(phi - phi_old), [], 'all') < tolerance
break;
end
end
Ex = -diff(phi, 1, 2) / dx;
Ey = -diff(phi, 1, 1) / dy;
[X, Y] = meshgrid(linspace(x_min, x_max, Nx), linspace(y_min, y_max, Ny));
figure;
contourf(X, Y, phi);
colorbar;
In this code, we define the dimensions of the rectangle and the number of grid points in each direction. We then initialize the potential matrix and set the boundary conditions. The Gauss-Seidel iteration is performed until convergence is achieved. Finally, we calculate and plot the electric field.
I hope this helps!

类别

在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by