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

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.

[VBBV]

It's recommended not to use standard MATLAB function names as variables. integral is standard MATLAB function

Answer 1

[SAI SRUJAN]

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!