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

Use finite difference method to solve for eigenvalue E

  1. Dec 7, 2012 #1
    Use finite difference method to solve for eigenvalue E from the following second order ODE:

    - y'' + (x2/4) y = E y

    I discretize the equation so that it becomes

    yi-1 - [2 + h2(x2i/4)] yi + yi+1 = - E h2 yi

    where xi = i*h, and h is the distance between any two adjacent mesh points.

    This is my Matlab code:

    Code (Text):
    clear all
    n = 27;
    h = 1/(n+1);
    voffdiag = ones(n-1,1);
    for i = 1:n
        xi(i) = i*h;
    mymat = -2*eye(n)-diag(((xi.^2).*(h^2)./4),0)+diag(voffdiag,1)+diag(voffdiag,-1);
    lam= -D/(h^2);
    fprintf(1,' The smallest eigenvalue is %g \n',lam(1));
    fprintf(1,' The second smallest eigenvalue is %g \n',lam(2));
    fprintf(1,' The 3rd eigenvalue is %g \n',lam(3));
    it returns

    Code (Text):
     The smallest eigenvalue is 9.92985
     The second smallest eigenvalue is 39.3932
     The 3rd eigenvalue is 88.0729
    Obviously, something wrong here, since the analytic solution should be
    E = N + 1/2 (for N = 0, 1, 2, 3...)
    The smallest eigenvalue should be 0.5, instead of 9.92985.

    I don't know whether my numerical solution agrees with the analytic solution or not, if I impose a boundary condition (ie. when x goes to infinity, y(x) should vanish to 0). And I don't know how to impose boundary condition. Please help, thank you very much!

    By the way, is there any another numerical way to find the eigenvalue E please?
    Last edited: Dec 7, 2012
  2. jcsd
  3. Dec 7, 2012 #2


    User Avatar
    Science Advisor
    Homework Helper

    Are you sure? Shouldn't the eigvenvalues depend on the differential equation and its boundary conditions, not on the number of data points you use to solve it numerically?

    Try a different value of n, and see if you get approximately the same eigenvalues. That is what I would expect, if your Matlab code is correct.

    Your first three eigenvalues from Matlab are almost in the ratio 1 : 4 : 9, which suggests they aren't "random numbers" but trying to tell you something.
  4. Dec 7, 2012 #3
    Thank you so much for your reply, AlephZero!

    Oh, I am sorry that I use the same letter to confuse people. In the code, n=27, which is the number of mesh points for discretized scheme; whereas E = n + 1/2 in which the n is just a finite interger starting from 0, in physics, this n is called quantum number. Since both n have different meanings, so I edit my main message by replacing E = n + 1/2 to E = N + 1/2.

    When I set n=177, my first three eigenvalues from Matlab are indeed approaching to the ratio 1 : 4 : 9. The following values are the first 6 eigenvalues.
    Code (Text):
    All the following values should times 10^5 in order to be correct value  
    Note that E = N+1/2 for N = 0, 1, 2, .... is indeed the analytic solution if we impose the boundary condition, that is, as x gos to infinity or negative infinity, y(x) should vanish to 0.

    However, the question is how to impose the boundary condition to my discretized scheme?
    Thank you so much!
  5. Dec 7, 2012 #4


    User Avatar
    Science Advisor
    Homework Helper

    I don't think you have told us what the boudary conditions are!

    I didn't read your code very carefully, but you seem to have something like y(0) = y(1) = 0.
  6. Dec 7, 2012 #5
    This problem looks a lot like "find the eigenfunctions of a quantum harmonic oscillator." (I had to write code for that assignment too.) If so, the usual boundary conditions are ##y(x) \to 0## as ##x \to \pm \infty##. You can't impose those boundary conditions on a finite-difference method because it has no way of representing infinitely-large values of ##x##.

    However, all is not lost. Your grid of ##x## values has to start and end somewhere. Let's call the least and greatest ##x## values ##x_L## and ##x_R##. The program ignores ##y(x)## if ##x## is "outside the grid." That's the same result you'd get if you used boundary conditions ##y(x) = 0## unless ##x_L \leq x \leq x_R##. If you choose ##x_L, x_R## far enough to the left and right, then this is a pretty good approximation to the true boundary conditions.

    My rule of thumb: the program adds a particle-in-a-box potential to whatever potential you assigned. So be sure to use a big box.
  7. Dec 8, 2012 #6
    I have told the boundary condition on my very first message:

    Hello, NegativeDept! Thank you very much for your reply! Although the program ignores y(x) if x is "outside the grid", this is not equivalent to the boundary condition. Because boundary condition says that as x grows larger and larger, y(x) should gradually become smaller and smaller, in other words, as x→∞, y(x) is still differentiable. However, the program ignoring y(x) if x is outside of grid means that y(x) can suddenly become 0 without any precursor, so that y(x) is not differentiable near the boundary.

    So imposing boundary condition is essential, and if we impose that, we will obtain a different result, very different from the result obtained by ignoring y(x) if x is outside of grid.

    We can choose the edge of the grid which is a finite and large value of x as pseudo-infinity, and impose pseudo-boundary condition. Still thank you, NegativeDept!

    Any another idea please?
  8. Dec 8, 2012 #7


    User Avatar
    Science Advisor
    Homework Helper

    Oops, I missed that.

    As well as just taking a big interval for x, you might try transforming the equation by letting v = 1/x say.

    Then ##0 \le v \le 1## corresponds to ##1 \le x \le \infty##. So you could piece together a solution for ##0 \ge v \ge -1##, ##-1 \le x \le 1##, ##1 \ge v \gt 0##.

    You would have to figure out how to set up the finite differences across the "boundary" between the x's and the v's.
  9. Dec 9, 2012 #8
    Oh, gosh! I just realize that physics forum does not support "copy and paste" so that the equations on my first message are totally messed up. Here is the modified version. I hope people can see and understand this:

    "Use finite difference method to solve for eigenvalue E from the following second order ODE:

    - y'' + ([itex]\frac{x^{2}}{4}[/itex]) y = E y

    I discretize the equation so that it becomes

    y[itex]_{i-1}[/itex] - [2 + h[itex]^{2}[/itex](x[itex]^{2}[/itex][itex]_{i}[/itex]/4)] y[itex]_{i}[/itex] + y[itex]_{i+1}[/itex] = - E h[itex]^{2}[/itex] y[itex]_{i}[/itex]

    where x[itex]_{i}[/itex] = i*h, and h is the distance between any two adjacent mesh points.

    The question is that I don't know how impose boundary condition which is as |x|[itex]\rightarrow[/itex][itex]\infty[/itex], y(x) [itex]\rightarrow[/itex]0."
  10. Dec 9, 2012 #9
    Thank you, Alephzero! Maybe your idea to change variable is useful and worth of try.

    However, this is just a common trick. I wish to see more detailed and specific solution.
    Last edited: Dec 9, 2012
  11. Dec 10, 2012 #10


    User Avatar
    Science Advisor
    Homework Helper

    Click on the link to "Rules" (which you agreed to when you signed up, even if you didn't read them!) at the top of the page, and see the section on "homework help guidelines".

    I gave you an outline of how to set up the problem this way. If you want a complete solution given to you, sorry, but PF is the wrong place to ask.

    If you try to solve it that way yourself, and get stuck, then of course you can ask some more specific questions about what is wrong.
  12. Dec 10, 2012 #11

    This is NOT a homework! If this is a homework, then this homework must be TOO easy! Also, homework usually requires students to use a specific numerical method to solve problem. However, here this question does not restrict people to come up with a different method to solve the problem. Any methods in addition to finite difference method, as long as it is numerical method, are perfectly fine!

    It is a personal interest to find a numerical way to solve elementary problem.

    I have already provided my attempted solution on my very first message, and got stuck on imposing the boundary value.
    Last edited: Dec 10, 2012
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook