# Find the electric field between 2 finite discs

• kirito
kirito
Homework Statement
there are two finite disc distance d from each other in the z axis potential difference between them V_zero find capacitance q and electric field between them
Relevant Equations
relaxation method , Q/C=V
I did make the problem simpler by looking at the the part from d/2 down the upper plate
here are my initial parameters I am making my size step be h since lowering it may make calculating harder
I am especially getting weird results for the field and capacitance

Code:
R = 0.1; % Radius of the plates in meters (10 cm)
d = 0.005; % Separation distance in meters (0.5 cm)
V0 = 1; % Potential difference in volts
h = d / 20; % Step size
rmax = R + 5 * d; % Maximum r
zmax = 10 * d; % Maximum z

% Create the grid
nr = ceil(rmax / h);
nz = ceil(zmax / h);
phi = zeros(nr, nz);

phi(1:R/h, d/(2*h)) = V0/2;% Set the potential on the positive plate
phi(end, :) =  0%; Boundary condition far from the plates
phi(:, 1) = 0 ;

% Relaxation method)
max_iterations = 10000;
tolerance = 1e-6;
for iter = 1:max_iterations

old_phi = phi;

for i = 2:nr-1
for j = 2:nz-1

if j == d/(2*h) && i < R/h
continue
end

phi(i, j) = (1/4) * (phi(i+1, j) * (1 + h / (2 * i)) + ...
phi(i-1, j) * (1 - h / (2 * i)) + ...
phi(i, j+1) + phi(i, j-1));
end
end
phi(:,1) = old_phi(:,2);
% Check for convergence
if max(max(abs(phi - old_phi))) < tolerance
fprintf('Converged after %d iterations\n', iter);
break;

end
end

% Plot the potential
figure;
imagesc((0:nr-1) * h, (0:nz-1) * h, flipud(phi'));
colorbar;
title('Electric Potential \phi');
xlabel('r (m)');
ylabel('z (m)');

% Compute the electric field
[B]% now here this seems like  I am calculating wrongly [/B]
[Er, Ez] = gradient(-phi, h);

% Compute the surface charge density at the edge of the plate
r_edge = ceil(R / h) + 1;
sigma = -Er(r_edge, :);

[B]% and the capacitance [/B]
% Note: We multiply by 2*pi*r to account for cylindrical coordinates
Q = 0;
for j = 1:nz-1
r = (j - 1) * h;
area = 2 * pi * r * h; % Area of the thin ring
Q = Q + sigma(j) * area;
end

% Compute the capacitance
% Compute the capacitance
C = Q / V0;
fprintf('Computed capacitance: %e F\n', C);

% Analytic capacitance for comparison
epsilon0 = 8.854187817e-12; % Vacuum permittivity
C_analytic = epsilon0 * pi * R^2 / d;
fprintf('Analytic capacitance: %e F\n', C_analytic);

% Plot the potential along z at r = 0
figure;
plot((0:nz-1) * h, phi(1, :);
title('Electric Potential \phi along z at r = 0');
xlabel('z (m)');
ylabel('\phi (V)');

% Plot the electric field
figure;
quiver((0:nr-1) * h, (0:nz-1) * h, Er', Ez');
title('Electric Field');
xlabel('r (m)');
ylabel('z (m)');

% Plot the surface charge density
figure;
plot((0:nz-1) * h, sigma);
title('Surface Charge Density \sigma at the Edge of the Plate');
xlabel('z (m)');
ylabel('\sigma (C/m^2)');

Last edited by a moderator:
kirito said:
phi(i+1, j) * (1 + h / (2 * i))
I am not seeing the logic behind this expression. Please explain it.
E.g., i is a dimensionless number but h has units of metres, so 1 + h / (2 * i) is dimensionally inconsistent.

haruspex said:
I am not seeing the logic behind this expression. Please explain it.
E.g., i is a dimensionless number but h has units of metres, so 1 + h / (2 * i) is dimensionally inconsistent.
you are correct I have preciously defined r=hi so if I have r_2=2h r_3=3h and so on then I said I will substitute it instead and forgot to delete the h from the upper term making it into 1+ 1/(2)

I did improve the code by a lot since I had a lot of confusion In application especially the relaxation method it is a bit hard for me visualise how it exactly works In a table I will add an updated version now

haruspex said:
I am not seeing the logic behind this expression. Please explain it.
E.g., i is a dimensionless number but h has units of metres, so 1 + h / (2 * i) is dimensionally inconsistent.
with the dimension inconsistency dealt with ,really shows how important it is to check units it took me until yesterday morning to notice it
an improved code with comments for a bit of clarity:
% Parameters
R = 0.1; % Radius of the plates in meters (10 cm)
d = 0.005; % Separation distance in meters (0.5 cm)
V0 = 1; % Potential difference in volts
h = d / 20; % Step size
rmax = R + 5 * d; % Maximum r
zmax = 10 * d; % Maximum z
epsilon0 = 8.854187817e-12; % Vacuum permittivity

% Create the grid
nr = ceil(rmax / h);
nz = ceil(zmax / h);
phi = zeros(nr, nz);

phi(1:ceil(R/h), ceil(d/(2*h))) = V0/2; % Set the potential on the positive plate
phi(end, :) =  0; % Boundary condition far from the plates
phi(:, 1) = 0;
phi(:, end) = 0;

% Relaxation method
max_iterations = 1000000;
tolerance = 1e-6;
for iter = 1:max_iterations
old_phi = phi;
for i = 2:nr-1
for j = 2:nz-1
if j == ceil(d/(2*h)) && i <= ceil(R/h)
continue
end
phi(i, j) = (1/4) * (old_phi(i+1, j) * (1 + 1 / (2 * i)) + ...
old_phi(i-1, j) * (1 - 1 / (2 * i)) + ...
old_phi(i, j+1) + old_phi(i, j-1));
end
end
phi(1,:) = phi(2,:);
% Check for convergence
if max(max(abs(phi - old_phi))) < tolerance
fprintf('Converged after %d iterations\n', iter);
break;
end
end

% Compute the electric field using central difference approximation
Er = zeros(nr, nz);
Ez = zeros(nr, nz);
for i = 2:nr-1
for j = 2:nz-1
Er(i, j) = -(phi(i+1, j) - phi(i-1, j)) / (2 * h);
Ez(i, j) = -(phi(i, j+1) - phi(i, j-1)) / (2 * h);
end
end

% Compute the surface charge density at the edge of the plate
%r_edge = ceil(R / h);
%sigma = epsilon0 * (Er(r_edge+1, :) - Er(r_edge-1, :)) / (2 * h);
% Compute the surface charge density at the edge of the plate
% Compute the surface charge density at the plate (at z = d/2)
z_plate = ceil(d / (2 * h)); % The index of z at the plate
sigma = epsilon0 * (Ez(:, z_plate+1) - Ez(:, z_plate-1)) / (2 * h);
% Compute the total charge on the plate by integrating sigma
% Note: We multiply by 2*pi*r to account for cylindrical coordinates
Q = 0;
for j = 1:nz-1
r = (j - 1) * h;
area = 2 * pi * r * h; % Area of the thin ring
Q = Q + sigma(j) * area;
end
fprintf('The total charge computed by integrating over the surface charged density:%e F\n', Q)
% Compute the capacitance
C = Q / V0;
fprintf('Computed capacitance: %e F\n', C);

% Analytic capacitance for comparison
C_analytic = epsilon0 * pi * R^2 / d;
fprintf('Analytic capacitance: %e F\n', C_analytic);

% Plot the potential
figure;
imagesc((0:nr-1) * h, (0:nz-1) * h, flipud(phi'));
colorbar;
title('Electric Potential \phi');
xlabel('r (m)');
ylabel('z (m)');

% Plot the potential along z at r = 0
figure;
plot((0:nz-1) * h, phi(1, :));
title('Electric Potential \phi along z at r = 0');
xlabel('z (m)');
ylabel('\phi (V)');

% Plot the electric field
figure;
quiver((0:nr-1) * h, (0:nz-1) * h, Er', Ez');
title('Electric Field');
xlabel('r (m)');
ylabel('z (m)');

% Plot the surface charge density
figure;
plot((0:nr-1) * h, sigma);
title('Surface Charge Density \sigma at the Edge of the Plate');
xlabel('r (m)');
ylabel('\sigma (C/m^2)');

% Numerical differentiation example
f = @cos;
x = 1;

d_f =@(x) -sin(x);
exact = d_f(x);
approx_1 = @(f, x, h) (f(x + h) - f(x)) ./ h;
approx_2 = @(f, x, h) (f(x + h) - f(x - h)) ./ (2 * h);
h = logspace(-1, -15, 100);

% Plot numerical differentiation errors
figure;
loglog(h, abs(exact - approx_1(f, x, h)), '*');
hold on;
loglog(h, abs(exact - approx_2(f, x, h)), '*');
hold off;
xlabel("log(h) from x=1");
ylabel("log(of the difference) from -sin(1)");
legend("normal approx", "symmetrical approximation");

This still looks wrong to me:
kirito said:
old_phi(i+1, j) * (1 + 1 / (2 * i)) + ... old_phi(i-1, j) * (1 - 1 / (2 * i)) + ... old_phi(i, j+1) + old_phi(i, j-1)
Why the 2*?

You seem to have fixed 0 potential at r=0 and r=rmax. Wouldn’t it be better to have it as linearly Vmin to Vmax at r=0, and maybe also at rmax?
Also, you set V0/2 at one plate, so shouldn’t you set -V0/2 at the other?

haruspex said:
This still looks wrong to me:

Why the 2*?

You seem to have fixed 0 potential at r=0 and r=rmax. Wouldn’t it be better to have those as linearly Vmin to Vmax? Or maybe that at r=0 and uniformly Vavg at rmax?
Also, you set V0/2 at one plate, so shouldn’t you set -V0/2 at the other?
the 2* because i was asked to take r as ih which is the step taken in the r direction
Why did not i set the potential on the other plate to -vo/2 because i only looked at the upper plate taking it only into account

kirito said:
the 2* because i was asked to take r as ih which is the step taken in the r direction
Why did not i set the potential on the other plate to -vo/2 because i only looked at the upper plate taking it only into account
Though i may i have applied it incorrectly but when i asked a teacher for a hint to simplify the calculations since the run time is long , he told me a hint is to only look at the upper plate and try to find symmetry for the problem in the r theta coordinates, so i looked at the upper plate and only calculated the potential in a quarter of the space assuming that from the symmetry for rotation the potential i would get in each of the 4 quarters is the same

kirito said:
Though i may i have applied it incorrectly but when i asked a teacher for a hint to simplify the calculations since the run time is long , he told me a hint is to only look at the upper plate and try to find symmetry for the problem in the r theta coordinates, so i looked at the upper plate and only calculated the potential in a quarter of the space assuming that from the symmetry for rotation the potential i would get in each of the 4 quarters is the same
haruspex said:
This still looks wrong to me:

Why the 2*?

You seem to have fixed 0 potential at r=0 and r=rmax. Wouldn’t it be better to have it as linearly Vmin to Vmax at r=0, and maybe also at rmax?
Also, you set V0/2 at one plate, so shouldn’t you set -V0/2 at the other?
I still have trouble understanding what potential to apply and for which condition, but i will reread this and try to understand how to apply what you are suggesting, i think it all comes down to the fact that though i know how to calculate relaxation I don’t understand how to apply Laplace conditions to most problems in question , i will be going to the office hours to see the code the lecturer wrote later on since the assignment time was due a while back , and i will read more about Laplace today

kirito said:
the 2* because i was asked to take r as ih which is the step taken in the r direction
For a given circle (r value and z value) you have four adjacent circles: (r+h,z), (r-h,z), (r,z-h), (r,z+h).
The contributions of these to the potential at (r,z) at the next time step must be weighted according their "volumes": r+h, r-h, r, r. Normalising, 1+h/r, 1-h/r, 1, 1. But you have 1+h/(2r) etc.
kirito said:
Why did not i set the potential on the other plate to -vo/2 because i only looked at the upper plate taking it only into account
Yes, sorry, I forgot you were modelling only the upper half.
But what about the potentials at r=0 and r=rmax? It is not right for those to be fixed at 0.

haruspex said:
For a given circle (r value and z value) you have four adjacent circles: (r+h,z), (r-h,z), (r,z-h), (r,z+h).
The contributions of these to the potential at (r,z) at the next time step must be weighted according their "volumes": r+h, r-h, r, r. Normalising, 1+h/r, 1-h/r, 1, 1. But you have 1+h/(2r) etc.

Yes, sorry, I forgot you were modelling only the upper half.
But what about the potentials at r=0 and r=rmax? It is not right for those to be fixed at 0.

The reason for the two because the formula i got when using symmetrical first derivative to approximate Laplace in spherical coordinates was this

kirito said:
View attachment 348616
The reason for the two because the formula i got when using symmetrical first derivative to approximate Laplace in spherical coordinates was this
For r=0 i thought its the condition so potential at z=0 is zero
For rmax i though i should take it as 0 since at infinite distance we consider potential zero

kirito said:
For r=0 i thought its the condition so potential at z=0 is zero
Surely at r=0 the potential should increase approximately linearly from 0 to V0/2.
kirito said:
For rmax i though i should take it as 0 since at infinite distance we consider potential zero
It depends whether rmax is sufficiently large in relation to the disc radius. You only have it 25% larger. Maybe experiment with that later when you have it basically working, i.e. see whether increasing rmax makes a significant difference.

kirito
haruspex said:
Surely at r=0 the potential should increase approximately linearly from 0 to V0/2.
Better would be to set each potential at r=0 at t+1 from the potential at r=h, same z, at time t.

haruspex said:
Better would be to set each potential at r=0 at t+1 from the potential at r=h, same z, at time t.
if you don't mind further explaining what this means , not sure if I get it , to set the potential at t+1 from r=h is bit unclear to me

kirito said:
if you don't mind further explaining what this means , not sure if I get it , to set the potential at t+1 from r=h is bit unclear to me
At each time step you set the potential at each (r,z) coordinate from its neighbours' potentials at the previous time step, but not for r=0. I'm saying set it at (0,z) from the value at (h,z) at the previous time step.
Or you could do an average involving (0,z+1), (0,z-1) too, but I don't think that changes much.

kirito
haruspex said:
At each time step you set the potential at each (r,z) coordinate from its neighbours' potentials at the previous time step, but not for r=0. I'm saying set it at (0,z) from the value at (h,z) at the previous time step.
Or you could do an average involving (0,z+1), (0,z-1) too, but I don't think that changes much.
your explanation is much appreciated, thank you

Replies
17
Views
911
Replies
64
Views
3K
Replies
1
Views
649
Replies
11
Views
465
Replies
6
Views
556
Replies
7
Views
564
Replies
1
Views
1K
Replies
6
Views
1K
Replies
5
Views
1K
Replies
2
Views
1K