I have code which calculates the electric potential due a single charge. The results are concentric circles (Correct results) but don't match the theoretical data

Here is the code clear all;

close all;

clc;

Nx = 11; %% Dimension of x axis

Ny = 11; %% Dimension of x axis

eps0 = 8.85e-12;

Q = 1; %% Charge Q (1C)

dx = 1; %% x step

dy = 1; %% y step

xq = ceil(Nx/2); %% x position of the charge

yq = ceil(Ny/2); %% y position of the charge

for i = 1 : Nx

    for j = 1 : Ny

        x = ( i * dx );

        y = ( j * dy );

        r = sqrt(( x-xq ).^2 + ( y-yq ).^2 );

        U(i,j) = Q ./ ( 4 .* pi .* eps0 .* r );

    end

end

% Contour Display for electric potential

figure

x = 1 : Nx;

y = 1 : Ny;

newpoints = 10;

contour( x/Nx, y/Ny, U, newpoints, 'linewidth', 3 );

% axis( [min(x/10) max(x/10) min(y/10) max(y/10)] );

colorbar( 'location', 'eastoutside', 'fontsize', 14 );

xlabel( 'r[m]', 'fontsize', 14 );

ylabel( 'U[V]', 'fontsize', 14 );

title( 'Electric Potential distribution in volts', 'fontsize', 14 );

h1 = gca;

set( h1, 'fontsize', 14);

fh1 = figure(1);

set( fh1, 'color', 'white' )

pbaspect( [1 1 1] );

%%% When I plot for instance the row 5th and column 5th to the end for the potential with respect to the distance 'r' and compare it with the theoretical plot of the formula (U = Q / 4 * pi * eps0 * r), the result doesn't match!!! I need your help guys! Thank you in advance.

Answer 1

[David Goodmanson]

Hi Matthew,

probably a bit late to be answering this now, but thanks to some editing done by John Kelly at Mathworks it provides an opportunity to delete an off-base comment that I made in 2018. I don't know what I was thinking.

The contour plot is correct, except that the axes are labeled incorrectly. The contour plot is

contour( x/Nx, y/Ny, U, newpoints, 'linewidth', 3 );

which means that the axes are the coordinates x/Nx and y/Ny, with U being indicated by the colors. The axes are not r and U.

As a detail I removed the Nx and Ny scaling with

contour(x,y, U, newpoints, 'linewidth', 3 );

so the axes are directly in meters,

xlabel( 'x[m]', 'fontsize', 14 );
ylabel( 'y[m]', 'fontsize', 14 );

There are nice circular isopotential lines in cartesian coordinates. Using datatips to pick a value off of the contour plot,

x = 4
y = 5
U = 4021 000 000

and since the charge is at (6,6) this is perfectly consistent with

U = Q/(4*pi*eps0*sqrt((4-6)^2+(5-6)^2))
   U = 4.0213e+09