# Need help getting this code to plot

3 views (last 30 days)
CAM on 2 Dec 2015
Answered: Ingrid on 2 Dec 2015
I want to plot the capacitance in this graph (cap) vs. the grid size to show how (in theory) it should get moe and more exact as grid size gets smaller. But I can't figure out how to do it. I tried putting it in the for loop, and outside the for loop, but have made no progress. Can someone help me figure out how to do this plot?
% --------------------------------------------------------------% Compute capacitance per unit length of
% a coaxial pair of rectangles
% -------------------------------------------------------------
function cap = capacitor(a, b, c, d, n, tol, rel)
% Arguments:
a = 0.015; %width of inner conductor
b = 0.015; %height of inner conductor
c = 0.02; %width of outer conductor
d = 0.02; %height of outer conductor
n = 10; %number of points in the x-direction (horizontal)
tol = 0.001; %relative tolerance for capacitance
rel = 1.7; %relaxation parameter
% (a good choice is 2-c/n, where c is about pi)
% Returns:
% cap = capacitance per unit length [pF/m]
% Make grids
h = 0.5*c/n; % Grid size
na = round(0.5*a/h); % Number of segments on ’a’
x = linspace(0,0.5*c,n+1); % Grid points along x-axis
m = round(0.5*d/h); % Number of segments on ’d’
mb = round(0.5*b/h); % Number of segments on ’b’
y = linspace(0,0.5*d,m+1); % Grid points along y-axis
% Initialize potential and mask array
f = zeros(n+1,m+1); % 2D-array with solution
mask = ones(n+1,m+1)*rel; % 2D-array with relaxation
% unchanged f(i,j)]
for i = 1:na+1;
for j = 1:mb+1;
f(i,j) = 1;
end
end
% Gauss Seidel iteration
oldcap = 0;
for iter = 1:1000000; % Maximum number of iterations
f = seidel(f,mask,n,m); % Perform Gauss-Seidel iteration
cap = gauss(n,m,h,f);% Compute the capacitance
if (abs(cap-oldcap)/cap<tol);
break % Stop if change in capacitance
% is sufficiently small
else
oldcap = cap; % Contiue until converged
end
end
str = sprintf('Number of iterations = %4i',iter); disp(str);
CAM on 2 Dec 2015
other relevant code:
function cap = gauss(n, m, h, f)
% Arguments:
% n = number of points in the x-direction (horizontal)
% m = number of points in the y-direction (vertical)
% h = cell size
% f = 2D-array with solution
% Returns:
% cap = capacitance per unit length [pF/m]
q=0;
for i = 1:n;
q=q+(f(i,m)+f(i+1,m))*0.5; % integrate along upper boundary
end
for j = 1:m;
q=q+(f(n,j)+f(n,j+1))*0.5; % integrate along right boundary
end
cap = q*4; % 4 quadrants
cap = cap*8.854187; % epsilon0*1e12 gives answer in pF/m
and
% --------------------------------------------------------------% Make one Seidel iteration
% -------------------------------------------------------------
function f = seidel(f, mask, n, m)
% Arguments:
% f = 2D-array with solution
% mask = 2D-array with relaxation
% n = number of points in the x-direction (horizontal)
% m = number of points in the y-direction (vertical)
% Returns:
% f = 2D-array with solution after Gauss-Seidel iteration
% Gauss seidel iteration
for i = 2:n
for j = 2:m
f(i,j) = f(i,j) + mask(i,j)* ...
(0.25*( f(i-1,j) + f(i+1,j) ...
+ f(i,j-1) + f(i,j+1)) - f(i,j));
end
end
% Symmetry on left boundary i-1 -> i+1
i=1;
for j = 2:m
f(i,j) = f(i,j) + mask(i,j)* ...
(0.25*( f(i+1,j) + f(i+1,j) ...
+ f(i,j-1) + f(i,j+1)) - f(i,j));
end
% Symmetry on lower boundary j-1 -> j+1
j=1;
for i = 2:n;
f(i,j) = f(i,j) + mask(i,j)* ...
(0.25*( f(i-1,j) + f(i+1,j) ...
+ f(i,j+1) + f(i,j+1)) - f(i,j));
end

Ingrid on 2 Dec 2015
it is not clear what you want to do but is it something like this?
function plotting(gridSizes)
cap = zeros(length(gridSizes),1);
for ii = 1:length(gridSizes)
cap(ii) = capacitor(gridSizes(ii));
end
plot(gridSizes,cap);
end
you could call this function with
gridSizes = 1:10;
plotting(gridSizes);