Help with Method of Moments, Surface Charge on Plate

Hello!

I am reading the book "The Method of Moments in Electromagnetics" by Walton C. Gibson, p. 39-43. (If anyone has it, I found a PDF online very easily.) Specifically, it is for a charged plate. Now, in the book he uses approximations for the integrals to obtain the surface charge. However, I want to see if I can do it by numerical calculation of the integrals (except the self terms). I get something close to his, namely, singularities near the edges and "dipping down" near the center. But it is only happening in one dimension! Here is my code (I have also attached the file so it is easier to read in Matlab, although I guess you can just copy and paste). If you run it, you will see what I mean. Hope this isn't too much.

clc
clear all
V = 1; %Voltage
L = 1; %Length of Plate in X and Y dimension
N = 10; %Number of segments to divide both X and Y directions into
n = 20; %Numer of divisions per segment in each X and Y direction for
      %integration over a patch. (So the whole plate is divided into an
      %NxN array of patches, and each patch is divided into an nxn array
      %for integration.)
Delta = L/N; %Length of patch on each side.
eo = 8.854e-12; %Permittivity of Free Space
v = V*4*pi*eo*ones(N^2,N^2); %Voltage vector multiplied by 4*pi*eo
%Calculate Source Coordinates (which are in the centers of the patches)
Xm = 0.5*(-L + Delta) : Delta : 0.5*(L-Delta);
Ym = 0.5*(-L + Delta) : Delta : 0.5*(L-Delta);
%Calculate integration limits for any patch (Observation coordinates)
for j = 1:N
  X{j} = Xm(j)-Delta/2 : Delta/n : Xm(j)+Delta/2;
  Y{j} = Ym(j)-Delta/2 : Delta/n : Ym(j)+Delta/2;
end
%Make an array of all of the XY coordinates of the source points
s = 1;
for k = 1:length(Xm)
  for l = 1:length(Ym)
      P(s,:) = [Xm(k),Ym(l)];
      s = s + 1;
  end
end
%Calculate the Z matrix components.
for M = 1:length(P)
  NN = 1;
  for J = 1:length(X)
      for K = 1:length(Y)
          for jj = 1:length(X{J})
              for kk = 1:length(Y{K})
                  z(jj,kk) = 1./(sqrt((P(M,1)-X{J}(jj)).^2+(P(M,2)-Y{K}(kk)).^2));
              end
          end
          Z(M,NN) = trapz(Y{K},trapz(X{J},z,1),2);
          NN = NN + 1;
      end
  end
end
%Replace all self terms with analytical evaluation of integral to avoid
%singularities.
for ZZ = 1 : length(Z)
  Z(ZZ,ZZ) = 8*Delta*asin(1);
end
%Calculate charge from Z matrix.
Q = v*Z^(-1)/(10^(-12));
%Get list of points such that the center point of each patch corresponds
%to the calculated charge on that patch.
l = 1;
for S = 1:N
  for s = 1:N
      POINTS(l,:) = [Xm(S),Ym(s),Q(S,s)];
      l = l + 1;
  end
end
scatter3(POINTS(:,1),POINTS(:,2),POINTS(:,3))
xlabel('X')
ylabel('Y')
zlabel('Z')

[Matthew Brandsema]

Note that for the section "Calculating the Z components" the way that works is, first I find the coordinates of 1 patch by dividing it into an n x n grid. This is done by using the jj and kk counting variables. I then integrate it (the capital Z term). I should have put more comments there.

Answer 1

[David Goodmanson]

Hi Matthew, I believe the basic problem has to do with the representation of voltage and charge. You have a 10x10 grid of points, so 100 elements. The Z matrix should be 100x100, which it is. Z is operating on a column vector containing all the (unknown) charge elements. That vector has to be related to the 10x10 square matrix for Q. This is usually done in reversible fashion by reading the elements in and out columnwise:

Q_one_column = Q(:)
Q = reshape(Q_one_column,10,10)   % back to square matrix, 10x10 in this case

It looks like that is what is happening with the for loops you have. Then you divide Z into the known v. However, v at the moment has 10000 values, not 100. If you go with

% voltage vector, same voltage on each
v = V*4*pi*eo*ones(N^2,1);   %Voltage vector multiplied by 4*pi*eo
Q = (Z\v)/1e-12;
Q = reshape(Q,N,N);
mesh(1:N,1:N,Q)        % quick, not correct x and y axes for demo purposess

you will get something that looks a lot better. (You are better qualified than I to judge if the scaling of Q is correct). The scatter3 plot also works at this point.

It's a bit tricky in this case because of the 10x10 ->100x100 ->10x10 aspect, but the meshgrid capability of Matlab can be really helpful as a good alternative to 'for' loops.

Most importantly, though, the code above does NOT use Z^(-1) or the equivalent inv(Z). You only want to find Q and are not making any explicit use of inv(Z). It is faster and more accurate to use the backslash operator, Q = Z \ v. I don't believe it's much of an exaggeration to say that the backslash operator was the genesis for Mathworks.

[Matthew Brandsema]

Makes sense! I really didn't expect to get any good answers as my question involves a bit of work to see where the error is (I would guess). Thanks so much!

Also, the values are indeed correct for Q.