Use finite difference method to solve for eigenvalue E

EigenCake

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:

clear all
n = 27;
h = 1/(n+1);
voffdiag = ones(n-1,1);
for i = 1:n
    xi(i) = i*h;
end
mymat = -2*eye(n)-diag(((xi.^2).*(h^2)./4),0)+diag(voffdiag,1)+diag(voffdiag,-1);
D=sort(eig(mymat),'descend');
lam= -D/(h^2);
spy(mymat)
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:

 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?

AlephZero

Quote by EigenCake

Obviously, something wrong here, since the analytic solution should be
E = n + 1/2 (for n = 0, 1, 2, 3...)

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.

EigenCake

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:

  
All the following values should times 10^5 in order to be correct value  
    0.0001
    0.0004
    0.0009
    0.0016
    0.0025
    0.0036

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!

AlephZero

Quote by EigenCake

However, the question is how to impose the boundary condition to my discretized scheme?

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.

NegativeDept

Quote by EigenCake

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!

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.

EigenCake

Quote by AlephZero

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

I have told the boundary condition on my very first message:

when x goes to infinity, y(x) should vanish to 0

---------------------------------------------------------------------------------------------------------

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 ≤x≤x R .

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.

You can't impose those boundary conditions on a finite-difference method because it has no way of representing infinitely-large values of x .

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?

AlephZero

Quote by EigenCake

I have told the boundary condition on my very first message:

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.

EigenCake

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."

EigenCake

Quote by AlephZero

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.

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.

AlephZero

Quote by EigenCake

I wish to see more detailed and specific solution.

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.

EigenCake

Quote by AlephZero

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".

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.

Similar discussions for: Use finite difference method to solve for eigenvalue E
Thread	Forum	Replies
Use finite difference method to solve for eigenvalue E in Matlab	Math & Science Software	0
Finite element method versus intergrated finite difference for complex geometries?	Differential Equations	6
Finite Difference method to solve PDEs	Differential Equations	1
Implicit Finite difference method to solve a Heat Transfer problem	Mechanical Engineering	0
stability conditions for finite difference method to solve coupled equations	Differential Equations	0