Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Nonuniform finite element method

  1. Feb 14, 2016 #1
    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:

    2ψi/∂x2 = 1/(Δx)2i-1i+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 = Δxi, 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?
     
  2. jcsd
  3. Feb 14, 2016 #2

    andrewkirk

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    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}$$
     
  4. Feb 15, 2016 #3
    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.
     
  5. Feb 15, 2016 #4

    andrewkirk

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    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}##.
     
  6. Feb 16, 2016 #5
    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);
     
  7. Feb 16, 2016 #6

    andrewkirk

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    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.
     
  8. Feb 17, 2016 #7
    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.
     
  9. Feb 17, 2016 #8

    andrewkirk

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    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:
    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.
     
  10. Feb 17, 2016 #9
    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?
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Nonuniform finite element method
  1. Finite element method (Replies: 1)

Loading...