-3
$\begingroup$

I know the charge density at the plates of a capacitor, from which the applied voltage and potential profile is to be calculated. Using the finite difference method in MATLAB script I solved the Poisson's equation that gives the correct result for all kind of continuous charge profiles. I checked constant (step function), linearly increasing/decreasing charge profile, in all cases numerical results agree with analytically calculated potential profile. But when I put the delta-function charge to study the original problem (in figure), it gives spooky results. (With $\rho=13.2$$\mu$C/$cm^2$, $\epsilon_r$=30, $x$=20nm, I would expect the voltage to be 9.92V,while I am getting 40nV!)

In the finite different approach:

$d^2V/dx^2$=-$\rho(x)$/$\epsilon$ in matrix form is $Au$=$h^2$$f$, where $f$=-[$\rho(i)$/$\epsilon$]' and $i$=1,2,......n are the discrete nodes. In my problem, I have $-\rho(1)$=$\rho(n)$=$13.2$$\mu$C/$cm^2$, and $\rho(i)=0$ elsewhere. Therefore I take $f$=$-\rho(i)$/$\epsilon$[-1 0 0 0 .......... 0 0 0 1]'. Then potential is obtained by $u$=$A^{-1}$$f$

My question is, what is the physical reason for such anomaly?

How can I solve this issue and make this code work? enter image description here enter image description here

Code:

clc;

clear all;

t=20;  %Thickness value in nm

epsil_r=30; 

epsil=epsil_r*8.854*10^-14; %In CGS unit


n=500;              %Number of grid points


x1=0; x2=t*10^-7; %Thickness in cm


x=linspace(x1,x2,n);h=x(2)-x(1);


rho=13.2*10^-6;     %Strength of delta-function charge, C/cm^2

A=diag(-2*ones(1,n))+diag(ones(1,n-1),1)+diag(ones(1,n-1),-1);

A(n,n)=-1;        %Neumann boundary at right node

f=-rho/(epsil)*[-1;zeros(n-2,1);1]; %Right hand side of Poisson.

u=inv(A)*[h^2*f];              %Potential
$\endgroup$
  • 1
    $\begingroup$ Would Computational Science be a better home for this question? $\endgroup$ – Qmechanic Apr 30 '17 at 16:40
  • $\begingroup$ I think people will be hesitant to discuss the physical behaviour of a simulation whose code they could not verify ... ;) $\endgroup$ – Sanya Apr 30 '17 at 18:59
  • $\begingroup$ "What is the physical reason for such anomaly?" Either your computer is broken, or (more likely) there is a bug in your code. If you want the code debugged without showing it to anybody you need a mind-reading forum, not SE. $\endgroup$ – alephzero Apr 30 '17 at 23:48
  • $\begingroup$ Dear altruist, I added the code in the question. I was wondering if there is anything fundamental for which delta-function charge can not be used, or it is all about my code. I am expecting your guideline. $\endgroup$ – Alam May 2 '17 at 9:32
  • $\begingroup$ I don't understand the code in this question. There's a bunch of definitions, but I don't see what actually acts on those objects. How you plot things is not important, what's important is how you solve things. $\endgroup$ – Emilio Pisanty May 2 '17 at 10:10
1
$\begingroup$

Your units don't add up. You're solving for the equation $$ \frac{\mathrm d^2 u}{\mathrm dx^2} = f $$ where $f = \rho/\epsilon$. If you take your solution at face value, then $u$ will have the dimensionality \begin{align} [u] & = [h^2 f] = [h^2\rho/\epsilon_0] = [Q/\epsilon_0]. \end{align} Since $[\epsilon_0]=[Q/(VL)]$, this means that $[u]=[VL]$, a voltage times length, which is obviously not what you intended. Since the units don't match up, it is meaningless to talk about inconsistencies in the numbers - if you did this in some other unit system you would get a different answer.

The problem is in your handling of the surface charge associated with the delta function. The right-hand side of the Poisson equation requires a legitimate volumetric charge density, or something dimensionally equivalent to that. If $\sigma$ is a surface charge density, with $[\sigma]=[Q/L^2]$, then $$\rho(x,y,z) = \sigma\delta(x-x_0)$$ is a correct volumetric charge density because $[\delta(x-x_0)] = [1/L]$, but your discretization procedure has dropped this fact out.

To fix this, I would suggest dropping the pretense that you have a true singular charge density and instead model it as a finite slab of thickness $d$ and volumetric charge density $\sigma/d$, where you can set $d$ to your discretization length if you want. Ultimately, the value of $d$ should be unimportant to your finite result - but it's crucial that you get your units right.

$\endgroup$
  • $\begingroup$ Thank you very much, its giving perfect result now. I just needed to convert the surface charge density into "per volume" unit. @Emilio Pisanty $\endgroup$ – Alam May 2 '17 at 16:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.