Nonuniform finite element method

[aaaa202]

I am solving some 2nd order differential equations using the finite element method. Doing so I represent the second order derivative at a given point as:

∂²ψ_i/∂x² = 1/(Δx)² (ψ_i-1+ψ_i+1-2ψ_i)

And solve the differential equation by setting up a matrix of N entries and solving for the eigenvectors. I guess you are all familiar with this approach. Now what I want to do is use a nonuniform mesh. To do so I have tried to simply let Δx = Δx_i, i.e. such that it may vary depending on the site i. But doing so I get some solutions that don't make sense, which I assume is because of the "discontinuity" in the mesh size.

How can I approach the problem of a nonuniform mesh size?

 

[andrewkirk]

The trouble is that using the factor $(\Delta x)^2$ assumes that $x_{j+1}-x_j=x_j-x_{j-1}$, where that is not in fact the case.

An appropriate formula that avoids that would be

$$\frac{\partial^2\psi}{\partial x^2}(x_j)=\frac{\frac{\psi_{j+1}-\psi_j}{x_{j+1}-x_{j}}-
\frac{\psi_{j}-\psi_{j-1}}{x_j-x_{j-1}}}
{(x_{j+1}-x_{j-1})/2}$$

 

[aaaa202]

hmm when I construct the matrix for the formula above I still get something whose eigenvectors diverge at the point where x_j changes discontinuely.

 

[andrewkirk]

Post the problem you are trying to solve, and the results you've obtained so far.
If you have code representing your calculations - eg if it's in MATLAB or R - post that too.

I don't know if this is relevant, because it should be a minor issue unless the mesh is very coarse, but the formula above is not exactly an approximation for the second derivative at $x_j$. It is an approximation for the second derivative at $\frac{x_{j+1}-x_{j-1}}{2}$.

 

[aaaa202]

Okay I have the code here. Basically I am simulating a metal coupled to a heterostructure. Since the fermi wavelength is much shorter in the metal I want to make a grid, which has a higher density of points here. To do so I write:

%%%%% Defining the grid
N=2000; %Number of grid points
N0=2; %Parameter that defines how long outside the actual well the grid should extend N0=1 only contains the well.
N_Al=200;
N_Semi=N-N_Al;
dx=N0*L/N_Semi;
dx_Al=L_Al/N_Al;
x=[dx_Al*(1:N_Al),dx*((N_Al+1):(N_Semi+N_Al))];

I then calculate the kinetic energy operator, which is the second derivative (times some factor, which we for now will take as constant). That means it is proportional to the follow matrix:

for i=2:N-1
K(i,i)=-2/(x(i+1)-x(i-1))*(1/(x(i+1)-x(i))+1/(x(i)-x(i-1)));
K(i+1,i)=2/(x(i+1)-x(i-1))*(1/(x(i+1)-x(i)));
K(i,i+1)=2/(x(i+1)-x(i-1))*(1/(x(i)-x(i-1)));
end
K(N,N)=K(N-1,N-1);
K(1,1)=K(2,2);
K(1,2)=K(2,3);
K(2,1)=K(2,3);

I think this should be the correct expression from the one you gave. The reason for the last four lines is to make it the correct matrix while keeping it of size NxN.
Now I add to this the potential V, which is a diagonal matrix defined by (j labels different iteration cycles, but that is a longer store):

V(x<L,j)=-aff_GaAs;
V(x<L_InAs+L_GaAs,j)=-aff_InAs;
V(x<L_GaAs)=-aff_GaAs;
V(x<L_Al)=mu-Ef_Al;

And finally solve the eigenvalue problem:

H=K+diag(V(:,j),0);
[psi,E0] = eig(H); % Psi is a matrix with eigenvectors and E is a matrix with eigenenergies in diagonal.
E=diag(E0); % Vector of Eigenvalues
groundstate(:,j)=psi(:,1);
firstexcited(:,j)=psi(:,2);

 

[andrewkirk]

I then calculate the kinetic energy operator, which is the second derivative (times some factor, which we for now will take as constant). That means it is proportional to the follow matrix:

for i=2:N-1
K(i,i)=-2/(x(i+1)-x(i-1))*(1/(x(i+1)-x(i))+1/(x(i)-x(i-1)));
K(i+1,i)=2/(x(i+1)-x(i-1))*(1/(x(i+1)-x(i)));
K(i,i+1)=2/(x(i+1)-x(i-1))*(1/(x(i)-x(i-1)));
end
K(N,N)=K(N-1,N-1);
K(1,1)=K(2,2);
K(1,2)=K(2,3);
K(2,1)=K(2,3);

I think this should be the correct expression from the one you gave.

My knowledge of MATLAB is only very basic, so perhaps I'm missing something.
But I can't see how this relates to what I wrote. The only input it appears to have is a vector of $x_i$ values. There is no sign of any function $\psi$ of which the second derivative is being taken.

 

[aaaa202]

Well that would be the eigenvector of the matrix I solve for. When the Schrödinger equation is discretized finding the eigenstates amounts to finding the eigenvectors of the matrix operator. In other words ψ is a vector of n entries with the ith entry being ψ_i.

 

[andrewkirk]

OK. That makes more sense out of it.
Looking at the code again in the light of that, I think there's an error in the following lines:

K(i+1,i)=2/(x(i+1)-x(i-1))*(1/(x(i+1)-x(i)));
K(i,i+1)=2/(x(i+1)-x(i-1))*(1/(x(i)-x(i-1)));

I think they should be

K(i,i-1)=2/(x(i+1)-x(i-1))*(1/(x(i)-x(i-1)));
K(i,i+1)=2/(x(i+1)-x(i-1))*(1/(x(i+1)-x(i)));

You'd need to consequently change your last line to

K(N,N-1)=K(N-1,N-2);

Try that and see how you go.

 

[aaaa202]

I think I figured out my problem. I had forgotten some constants that made the problems dimensions wrong. Now suppose I want to add the fact that the effective mass in the kinetic energy operator should also be position dependent. Then my differential operator would be:
d/dx(1/m*(x) d/dx)

In discretized form how would this look when I compare to your expression for the second derivative. Would I just have 1/m*(i+1) and 1/m*(i-1) in front of the two terms?