Explanation of parallel plates capacitor implementation with finite element methods and with in-homogenous domain.

function cap

close all;

hx = 0.002;

vx = 0:hx:0.1;

hy = 0.002;

h = hx;

vy = 0:hy:0.1;

nx = length(vx);

ny = length(vy);

n = nx*ny;

A = zeros(n);

b = zeros(n,1);

eps1 = 4

eps2 = 14

for ix = 1:nx

for jy = 1:ny

i = (jy-1) * nx + ix; % the global index of the ix,jy vertex

x = vx(ix); % the geometrical x coordinate

y = vy(jy); % the geometrical y coordinate

if ix == 1

% (1)

A(i,i) = 1;

b(i) = 0;

elseif ix == nx

% (2)

A(i,i) = 1;

b(i) = 10;

elseif jy == 1

if ix < (nx+1)/2

% (3)

A(i,i) = eps1*(-(2/hx^2)-(2/hy^2));

A(i,i+1) = eps1*((1/hx^2)-(1/2*x*hx));

A(i,i-1) = eps1*((1/hx^2)-(1/2*x*hx));

A(i,i+nx) = eps1*(1/hy^2);

elseif ix > (nx+1)/2

% (4)

A(i,i) = eps2*(-(2/hx^2)-(2/hy^2));

A(i,i+1) = eps2*((1/hx^2)-(1/2*x*hx));

A(i,i-1) = eps2*((1/hx^2)-(1/2*x*hx));

A(i,i+nx) = eps2*(1/hy^2);

else

% (5)

A(i,i) = (eps1+eps2)*(-(2/hx^2)-(2/hy^2));

A(i,i+1) = eps2*((1/hx^2)-(1/2*x*hx));

A(i,i-1) = eps1*((1/hx^2)-(1/2*x*hx));

A(i,i+nx) = (eps1+eps2)*(1/hy^2);

end

elseif jy == ny

% (6)

A(i,i) = eps2*(-(2/hx^2)-(2/hy^2));

A(i,i+1) = eps2*((1/hx^2)-(1/2*x*hx));

A(i,i-1) = eps2*((1/hx^2)-(1/2*x*hx));

A(i,i-nx) = eps2*(2/hy^2);

elseif jy < (ny+1)/2 && ix < (nx+1)/2

% (7)

A(i,i) = eps1*(-(2/hx^2)-(2/hy^2));

A(i,i+1) = eps1*((1/hx^2)-(1/2*x*hx));

A(i,i-1) = eps1*((1/hx^2)-(1/2*x*hx));

A(i,i-nx) = eps1*(1/hy^2);

A(i,i+nx) = eps1*(1/hy^2);

elseif jy > (ny+1)/2 || ix > (nx+1)/2

% (8)

A(i,i) = eps2*(-(2/hx^2)-(2/hy^2));

A(i,i+1) = eps2*((1/hx^2)-(1/2*x*hx));

A(i,i-1) = eps2*((1/hx^2)-(1/2*x*hx));

A(i,i+nx) = eps2*(1/hy^2);

A(i,i-nx) = eps2*(1/hy^2);

elseif jy < (ny+1)/2

% (9)

A(i,i) = (eps1+eps2)*(-(2/hx^2)-(2/hy^2));

A(i,i+1) = eps2*((1/hx^2)-(1/2*x*hx));

A(i,i-1) = eps1*((1/hx^2)-(1/2*x*hx));

A(i,i+nx) = ((eps1+eps2)/2)*(1/hy^2);

A(i,i-nx) = ((eps1+eps2)/2)*(1/hy^2);

elseif ix < (nx+1)/2

% (10)

A(i,i) = (eps1+eps2)*(-(2/hx^2)-(2/hy^2));

A(i,i+1) = ((eps1+eps2)/2)*((1/hx^2)-(1/2*x*hx));

A(i,i-1) = ((eps1+eps2)/2)*((1/hx^2)-(1/2*x*hx));

A(i,i+nx) = eps2*(1/hy^2);

A(i,i-nx) = eps1*(1/hy^2);

else

% (11)

A(i,i) = (eps1+(eps2*3))*(-(2/hx^2)-(2/hy^2));

A(i,i+1) = eps2*((1/hx^2)-(1/2*x*hx));

A(i,i-1) = ((eps1+eps2)/2)*((1/hx^2)-(1/2*x*hx));

A(i,i+nx) = eps2*(1/hy^2);

A(i,i-nx) = ((eps1+eps2)/2)*(1/hy^2);

end

end

end

% solve the problem

u = A\b;

s = reshape(u,nx,ny);

[X,Y] = meshgrid(vx,vy);

surf(X,Y,s)

end

[Shubham Khatri]

Hello, Can you please elaborate on what exactly you are asking?

Answer 1

[Shubham Khatri]

Hello,

As far as the code is concerned, Please remove the last 'end' from the code for it to work as in the code below.

close all;
    hx = 0.002;
    vx = 0:hx:0.1; 
    hy = 0.002;
    h = hx;
    vy = 0:hy:0.1; 
    nx = length(vx);
    ny = length(vy);
    n = nx*ny;
    A = zeros(n);
    b = zeros(n,1);
    
    eps1 = 4
    eps2 = 14    
    for ix = 1:nx
        for jy = 1:ny
            i = (jy-1) * nx + ix; % the global index of the ix,jy vertex
            x = vx(ix); % the geometrical x coordinate
            y = vy(jy); % the geometrical y coordinate
            
            if ix == 1
               % (1)
               A(i,i) = 1;
               b(i) = 0;
            elseif ix == nx
               % (2)
               A(i,i) = 1;
               b(i) = 10;
            elseif jy == 1
                    if ix < (nx+1)/2
                        % (3)  
                       
                        A(i,i) = eps1*(-(2/hx^2)-(2/hy^2));
                        A(i,i+1) = eps1*((1/hx^2)-(1/2*x*hx));
                        A(i,i-1) = eps1*((1/hx^2)-(1/2*x*hx));
                        A(i,i+nx) = eps1*(1/hy^2);
                    elseif ix > (nx+1)/2
                        % (4)
                        A(i,i) = eps2*(-(2/hx^2)-(2/hy^2));
                        A(i,i+1) = eps2*((1/hx^2)-(1/2*x*hx));
                        A(i,i-1) = eps2*((1/hx^2)-(1/2*x*hx));
                        A(i,i+nx) = eps2*(1/hy^2);                            
                    else
                        % (5)
                        A(i,i) = (eps1+eps2)*(-(2/hx^2)-(2/hy^2));
                        A(i,i+1) = eps2*((1/hx^2)-(1/2*x*hx));
                        A(i,i-1) = eps1*((1/hx^2)-(1/2*x*hx));
                        A(i,i+nx) = (eps1+eps2)*(1/hy^2); 
                    end
            elseif jy == ny
                % (6)
                A(i,i) = eps2*(-(2/hx^2)-(2/hy^2));
                A(i,i+1) = eps2*((1/hx^2)-(1/2*x*hx));
                A(i,i-1) = eps2*((1/hx^2)-(1/2*x*hx));
                A(i,i-nx) = eps2*(2/hy^2);                            
            elseif jy < (ny+1)/2 && ix < (nx+1)/2
                    % (7)
                    A(i,i) = eps1*(-(2/hx^2)-(2/hy^2));
                    A(i,i+1) = eps1*((1/hx^2)-(1/2*x*hx));
                    A(i,i-1) = eps1*((1/hx^2)-(1/2*x*hx));
                    A(i,i-nx) = eps1*(1/hy^2);
                    A(i,i+nx) = eps1*(1/hy^2);
                elseif jy > (ny+1)/2 || ix > (nx+1)/2
                    % (8)
                    A(i,i) = eps2*(-(2/hx^2)-(2/hy^2));
                    A(i,i+1) = eps2*((1/hx^2)-(1/2*x*hx));
                    A(i,i-1) = eps2*((1/hx^2)-(1/2*x*hx));
                    A(i,i+nx) = eps2*(1/hy^2);                            
                    A(i,i-nx) = eps2*(1/hy^2); 
                elseif jy < (ny+1)/2 
                    % (9)
                    A(i,i) = (eps1+eps2)*(-(2/hx^2)-(2/hy^2));
                    A(i,i+1) = eps2*((1/hx^2)-(1/2*x*hx));
                    A(i,i-1) = eps1*((1/hx^2)-(1/2*x*hx));
                    A(i,i+nx) = ((eps1+eps2)/2)*(1/hy^2);                            
                    A(i,i-nx) = ((eps1+eps2)/2)*(1/hy^2); 
                elseif ix < (nx+1)/2  
                    % (10)
                    A(i,i) = (eps1+eps2)*(-(2/hx^2)-(2/hy^2));
                    A(i,i+1) = ((eps1+eps2)/2)*((1/hx^2)-(1/2*x*hx));
                    A(i,i-1) = ((eps1+eps2)/2)*((1/hx^2)-(1/2*x*hx));
                    A(i,i+nx) = eps2*(1/hy^2);                            
                    A(i,i-nx) = eps1*(1/hy^2); 
            else
                    % (11)
                    A(i,i) = (eps1+(eps2*3))*(-(2/hx^2)-(2/hy^2));
                    A(i,i+1) = eps2*((1/hx^2)-(1/2*x*hx));
                    A(i,i-1) = ((eps1+eps2)/2)*((1/hx^2)-(1/2*x*hx));
                    A(i,i+nx) = eps2*(1/hy^2);                            
                    A(i,i-nx) = ((eps1+eps2)/2)*(1/hy^2); 
                end
            end
        end
    
    
    % solve the problem
    u = A\b;
    s = reshape(u,nx,ny);
    [X,Y] = meshgrid(vx,vy);
    surf(X,Y,s)

Hope it helps