the use of field graphing tools

K = 8.99e9;

% Input the number of charges
n = input('Enter number of charges: ');

Error using input
Support for user input is required, which is not available on this platform.

% Initialize arrays to store charge values and coordinates
q = zeros(1, n);
x = zeros(1, n);
y = zeros(1, n);

minx = Inf;
maxx = -Inf;
miny = Inf;
maxy = -Inf;
% Input charge values
for i = 1:n
    q(i) = input(['Enter the charge in coulombs for charge ', num2str(i), ': ']);
end

% Input coordinates of the charges
for i = 1:n
    x(i) = input(['Enter x coordinate for charge ', num2str(i), ': ']);
    y(i) = input(['Enter y coordinate for charge ', num2str(i), ': ']);
    if x(i) > maxx
        maxx = x(i);
    end
    if x(i) < minx
        minx = x(i);
    end
    if y(i) > maxy
        maxy = y(i);
    end
    if y(i) < miny
        miny = y(i);
    end
end

% Define the grid for visualization
[X, Y] = meshgrid((minx-1):0.1:(maxx+1), (miny-1):0.1:(maxy+1));
Ex = zeros(size(X));
Ey = zeros(size(X));



% i tried to improve the out put [X, Y] = meshgrid((minx-1):0.01:(maxx+1), (miny-1):0.01:(maxy+1));
%Ex = zeros(size(X));
%Ey = zeros(size(X));
% yet it also lead to the unablitiy to see clearly if there is for example
% a way to change the arrows size or ask the computer for optimal view
% since i am viewing all possible options and most end up badly

% Calculate the electric field components
for i = 1:n
    R = ((X - x(i)).^2 + (Y - y(i)).^2 ).^1.5; % Distance cubed
    Ex = Ex + K * q(i) * (X - x(i)) ./ R;
    Ey = Ey + K * q(i) * (Y - y(i)) ./ R;
end

% Plot the electric field vectors
figure;
quiver(X, Y, Ex, Ey);
xlabel('x (m)');
ylabel('y (m)');
title('Electric Field Vectors');
axis equal;
grid on;

% Plot the potential distribution (if needed)
% Potential at each point
V = zeros(size(X));
for i = 1:n
    R = sqrt((X - x(i)).^2 + (Y - y(i)).^2 ); % Distance
    V = V + K * q(i) ./ R;
end

figure;
contourf(X, Y, V, 50, 'LineColor', 'none');
colorbar;
xlabel('x (m)');
ylabel('y (m)');
title('Potential Distribution');
axis equal;
grid on;

i wrote the following code and i am trying to improve a couple of things

1 i want the graph to fit the boundaries so it wont draw a grid that has field lines in parts and the rest is grid

2 when i enter a negative point charge it gets problamatic this is physically correct but i want to try and write something to deal with it

3 the field lines are write but they do not really reflect the behavoir well especially when i write something like 3 + and -1 on a destince 1 from each other they look just like 1 and -1 a bit better in the voltage aspect

now this is for improvment that seem okayish yet from my trial and error i could not implement

now the problamtic request is for the next part of the code i am trying to write i am trying to draw a visullization of the field between 2 finite planes

i started this part by writing an approximation for the derivative and checking that it works

h = logspace(-1,-15,100);
real_d = -sin(1);
real_d2 = -sin(1+h);
di_d = (cos(1+h)-cos(1))./h ;
simt_d = (cos(1+h)-cos(1-h))./(2*h) ;
plot (log10(h) , log10(real_d - di_d),'*')
hold on
plot (log10(h) , log10(real_d - simt_d),'*')
hold o
xlabel ("log(h) from x=1")
ylabel ("log(di) from -sin(1)")

i derived the equations i need

yet i am having a hard time trying to implment them

% Constants
d = 0.5; % distance between plates in cm
R = 10.0; % radius of the plates in cm
V = 1.0; % potential difference in V
h = d / 20; % step size in cm
r_max = R + 5 * d; % maximum r value in cm
z_max = 10 * d; % maximum z value in cm
tolerance = 0.0001; % convergence criterion (0.01%)

% Define the grid
r_points = floor(r_max / h) + 1;
z_points = 2 * floor(z_max / h) + 1;
phi = zeros(r_points, z_points);

% Set boundary conditions for the finite plates
for i = 1:floor(R / h) + 1
    % Positive plate at z = d/2
    phi(i, floor(z_max / h) + 1 + floor(d / (2 * h))) = 0.5 * V;
    % Negative plate at z = -d/2
    phi(i, floor(z_max / h) + 1 - floor(d / (2 * h))) = -0.5 * V;
end

% Relaxation method
max_diff = tolerance + 1; % Initialize max difference
while max_diff > tolerance
    max_diff = 0;
    phi_new = phi;
    for i = 2:r_points-1
        for j = 2:z_points-1
            % Skip the points on the plates
            if (i * h <= R) && (j == floor(z_max / h) + 1 + floor(d / (2 * h)) || j == floor(z_max / h) + 1 - floor(d / (2 * h)))
                continue; 
            end
            phi_new(i, j) = 0.25 * ( ...
                phi(i + 1, j) * (1 + h / (2 * (i * h))) + ...
                phi(i - 1, j) * (1 - h / (2 * (i * h))) + ...
                phi(i, j + 1) + ...
                phi(i, j - 1));
            max_diff = max(max_diff, abs(phi_new(i, j) - phi(i, j)));
        end
    end
    phi = phi_new;
end

% Plotting the potential
r = 0:h:r_max;
z = -z_max:h:z_max;
[R, Z] = meshgrid(r, z);
figure;
contourf(R, Z, phi', 50, 'LineColor', 'none');
colorbar;
xlabel('r (cm)');
ylabel('z (cm)');
title('Electric Potential in the r-z Plane');

the drawing seems weird too if some one can give a bit of guidenss on the use of relaxation

[dareen]

i am trying to implement something along thoese lines but i have discs instead of plates , to reach this objective i started from easier problems like superposition and tried to get used to drawing field lines not there yet , and here it seemed like a bit of a jump to try to do this right after but i would like to learn

[dareen]

i also skiped points on the disc planes because of discontinuty and such

[Walter Roberson]

(Executing the second part of your code)

% Constants
d = 0.5; % distance between plates in cm
R = 10.0; % radius of the plates in cm
V = 1.0; % potential difference in V
h = d / 20; % step size in cm
r_max = R + 5 * d; % maximum r value in cm
z_max = 10 * d; % maximum z value in cm
tolerance = 0.0001; % convergence criterion (0.01%)

% Define the grid
r_points = floor(r_max / h) + 1;
z_points = 2 * floor(z_max / h) + 1;
phi = zeros(r_points, z_points);

% Set boundary conditions for the finite plates
for i = 1:floor(R / h) + 1
    % Positive plate at z = d/2
    phi(i, floor(z_max / h) + 1 + floor(d / (2 * h))) = 0.5 * V;
    % Negative plate at z = -d/2
    phi(i, floor(z_max / h) + 1 - floor(d / (2 * h))) = -0.5 * V;
end

% Relaxation method
max_diff = tolerance + 1; % Initialize max difference
while max_diff > tolerance
    max_diff = 0;
    phi_new = phi;
    for i = 2:r_points-1
        for j = 2:z_points-1
            % Skip the points on the plates
            if (i * h <= R) && (j == floor(z_max / h) + 1 + floor(d / (2 * h)) || j == floor(z_max / h) + 1 - floor(d / (2 * h)))
                continue; 
            end
            phi_new(i, j) = 0.25 * ( ...
                phi(i + 1, j) * (1 + h / (2 * (i * h))) + ...
                phi(i - 1, j) * (1 - h / (2 * (i * h))) + ...
                phi(i, j + 1) + ...
                phi(i, j - 1));
            max_diff = max(max_diff, abs(phi_new(i, j) - phi(i, j)));
        end
    end
    phi = phi_new;
end

% Plotting the potential
r = 0:h:r_max;
z = -z_max:h:z_max;
[R, Z] = meshgrid(r, z);
figure;
contourf(R, Z, phi', 50, 'LineColor', 'none');
colorbar;
xlabel('r (cm)');
ylabel('z (cm)');
title('Electric Potential in the r-z Plane');

[dareen]

i have an idea for impletaion i set the boundries down as zero and zero in the first column i sart calculating from a guess then updating the first column with the new calculated potintial and so on yet did not reach satisfactory results

[dareen]

@Walter Roberson yeah i am mostly not closed on if i am implementing this part correctly

phi_new(i, j) = 0.25 * ( ...
    phi(i + 1, j) * (1 + h / (2 * (i * h))) + ...
    phi(i - 1, j) * (1 - h / (2 * (i * h))) + ...
    phi(i, j + 1) + ...
    phi(i, j - 1));
max_diff = max(max_diff, abs(phi_new(i, j) - phi(i, j)));
end
end
phi = phi_new; 

i derivied the expression then checked it online tried to write a code that updates the potential each run until it starts to stablize yet still trying to get the hang of the implementaion , i even tried to draw field line to feel better about how this looks yet the field lines were compact and hard to see anything thats happening same thing with that additional part one i did to get used to this is a new field in matlan, the results felt like i only had 2 plates and line between them looking now it seems density out side is diffrent which is a good sign , but with harder to visulize fields it is harder to feel confident with the lines i get when in part one they did not really give the results i amied for excactly for things like e3 and -e ,

Answer 1

[Rahul]

Hi Dareen,

I understand that you are working with modeling point charges to more advanced setups like simulating the field between two finite charged disks, where you are encountering a few challenges regarding visualization accuracy, boundary behavior, and grid setup

The issue with the point charge simulation could be with cluttered grids and exaggerated field vectors near charges, especially negatives. To optimize this, you can consider tightening the plotting range based on the charge locations. For example, adding a small padding like 0.5 around the extrema could be helpful, as shown below:

[X, Y] = meshgrid((minx - 0.5):0.05:(maxx + 0.5), (miny - 0.5):0.05:(maxy + 0.5));

This can help center the field. In case the field vectors near charges are cluttering the plot, you can normalize the field for better directionality as shown below:

magnitude = sqrt(Ex.^2 + Ey.^2);
Ex_unit = Ex ./ magnitude;
Ey_unit = Ey ./ magnitude;
quiver(X, Y, Ex_unit, Ey_unit, 
0.5);  % Adjust the 0.5 scale as needed

Alternatively, you can use 'streamslice(X, Y, Ex, Ey)' which provides a better visual of the field lines and tends to behave better when many charges interact.

The relaxation phase includes a 2D cylindrical grid to solve Laplace’s equation between two disk electrodes. We can skip over the electrode boundaries (where discontinuity is a problem), and a modified stencil can be used to account for cylindrical symmetry. Although, the following approach for generating the mesh grid for plotting suggests the axes might be flipped:

r = 0:h:r_max;
z = -z_max:h:z_max;
[R, Z] = 
meshgrid(r, z);

To resolve this, you can switch to the following approach which eliminates the need for transposing 'phi':

[Z, R] = meshgrid(z, r);
contourf(Z, R, phi, 
50, 'LineColor', 'none');

You can also improve the look by adding:

shading interp;
colormap('jet');
axis equal;

Once the charge potential is stable, the electric field components can be visualized using the 'gradient' function as shown below:

[Er, Ez] = gradient(-phi, h);
quiver(Z, R, Ez, Er);  % Note the axis swap again

This can help in getting a result closer to the curved field lines seen in the provided image. If the iteration is slow, you can consider using Successive Over-Relaxation (SOR) to accelerate convergence, especially on denser grids.

To know more about the usage of various functions used in the given code snippets, you can refer to the following documentation links:

• 2-D and 3-D grids using 'meshgrid' function: https://www.mathworks.com/help/matlab/ref/meshgrid.html
• Numerical gradient using 'gradient' function: https://www.mathworks.com/help/matlab/ref/gradient.html
• Quiver plot using 'quiver' function: https://www.mathworks.com/help/matlab/ref/quiver.html

Hope this helps!