Use finite difference method to solve for eigenvalue E


by EigenCake
Tags: difference, eigenvalue, finite, method, solve
EigenCake
EigenCake is offline
#1
Dec7-12, 02:33 PM
P: 13
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:

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

 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?
Phys.Org News Partner Science news on Phys.org
Better thermal-imaging lens from waste sulfur
Hackathon team's GoogolPlex gives Siri extra powers
Bright points in Sun's atmosphere mark patterns deep in its interior
AlephZero
AlephZero is offline
#2
Dec7-12, 06:44 PM
Engineering
Sci Advisor
HW Helper
Thanks
P: 6,344
Quote Quote by EigenCake View Post
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
EigenCake is offline
#3
Dec7-12, 08:02 PM
P: 13
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.
  
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
AlephZero is offline
#4
Dec7-12, 11:02 PM
Engineering
Sci Advisor
HW Helper
Thanks
P: 6,344

Use finite difference method to solve for eigenvalue E


Quote Quote by EigenCake View Post
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
NegativeDept is offline
#5
Dec7-12, 11:31 PM
P: 135
Quote Quote by EigenCake View Post
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
EigenCake is offline
#6
Dec8-12, 10:23 PM
P: 13
Quote Quote by AlephZero View Post
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
AlephZero is offline
#7
Dec8-12, 11:12 PM
Engineering
Sci Advisor
HW Helper
Thanks
P: 6,344
Quote Quote by EigenCake View Post
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
EigenCake is offline
#8
Dec9-12, 10:39 PM
P: 13
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
EigenCake is offline
#9
Dec9-12, 10:42 PM
P: 13
Quote Quote by AlephZero View Post
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
AlephZero is offline
#10
Dec10-12, 01:20 PM
Engineering
Sci Advisor
HW Helper
Thanks
P: 6,344
Quote Quote by EigenCake View Post
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
EigenCake is offline
#11
Dec10-12, 02:41 PM
P: 13
Quote Quote by AlephZero View Post
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.


Register to reply

Related Discussions
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