# [MATLAB] Searching for specific eigenvalues with Newton iteration

• Jan 24th 2009, 05:14 AM
IgorMele
[MATLAB] Searching for specific eigenvalues with Newton iteration
Hi!
I would really appreciate if somebody could help me with my problem. I
have to calculate n-th eigenvalue of equation:

y'' = lambda * y

First I made a system of 4 equations:
http://www.shrani.si/t/1I/wQ/2shdGQhX/eigenvalues.jpg

then I found solution with RungeKutta method:
http://www.shrani.si/t/3W/Wh/3ywLKjI8/matlab.jpg

I have problems with Newton iteration. What to write in matlab.

Function:
Code:

f = inline('[-y(2); lam*y(1); -y(4); lam*y(3)+y(1)]','x','y','lam');
Here is my RK4:
Code:

function [x1, y1, spr, values] = RungeKutta4(f, x0, y0, xb, lam, N) %[x1, y1, lam, spr, values] = RungeKutta4(f, x0, y0, xb, lam, N) %y0 =[0 1 0 0]' %lam = lambda %N j=number of steps for RK4 %spr = number of changes: spr(++ ----++ - ++) = 4 values = []; values(1)=y0(1); h=(xb-x0)/N; spr = 0; for i = 1:N     k1 = h*feval(f, x0, y0, lam);     k2 = h*feval(f, x0+(h/2), y0+(k1/2),lam);     k3 = h*feval(f, x0+(h/2), y0+(k2/2),lam);     k4 = h*feval(f, x0+h, y0+k3,lam);     a=y0(1);     y0 = y0 + (k1 + 2*k2 + 2*k3 + k4)/6;     x0 = x0 + h;     values(i+1)=y0(1);     if sgn(a)~=sgn(y0(1))         spr = spr +1;     end     %hold on     %plot(x0,y0(1),'-x')     %plot(x0,y0(3),'-*') end
I made also program for searching interval where n-th eigenvalue is:
Code:

function root = BISECTION(f,xa,y0,xb,lam,N,c,d,whichLambda) root=(c+d)/2; while (root~=c) & (root~=d)       if whichLambda>SIGN(f,xa,root,lam,N)             c=root;     else         d=root;     end     root=(c+d)/2; end
SIGN counts number of sign changes in interval [xa, root]

And here is problem...newton iteration. My program:

Code:

function lam = fN(f, x0, y0, xb, lam, N, n) %lam = fN(f, x0, y0, xb, lam, N, n) %n = number of steps for Newton's iteration %N = number of steps for RK4 lam1=16.1; for i = 1:n     [x1, y1, spr, vr] = RungeKutta4(f, x0, y0, xb, lam, N);     lam1 = lam1 - (y1(1))/(y1(3)); end
I'm not getting any good results for eigenvalue (lam1).

Thanks for any help...

Igor
• Jan 24th 2009, 08:04 AM
Constatine11
Quote:

Originally Posted by IgorMele
Hi!
I would really appreciate if somebody could help me with my problem. I
have to calculate n-th eigenvalue of equation:

y'' = lambda * y

First I made a system of 4 equations:
http://www.shrani.si/t/1I/wQ/2shdGQhX/eigenvalues.jpg

then I found solution with RungeKutta method:
http://www.shrani.si/t/3W/Wh/3ywLKjI8/matlab.jpg

I have problems with Newton iteration. What to write in matlab.

Function:
Code:

f = inline('[-y(2); lam*y(1); -y(4); lam*y(3)+y(1)]','x','y','lam');
Here is my RK4:
Code:

function [x1, y1, spr, values] = RungeKutta4(f, x0, y0, xb, lam, N) %[x1, y1, lam, spr, values] = RungeKutta4(f, x0, y0, xb, lam, N) %y0 =[0 1 0 0]' %lam = lambda %N j=number of steps for RK4 %spr = number of changes: spr(++ ----++ - ++) = 4 values = []; values(1)=y0(1); h=(xb-x0)/N; spr = 0; for i = 1:N       k1 = h*feval(f, x0, y0, lam);     k2 = h*feval(f, x0+(h/2), y0+(k1/2),lam);     k3 = h*feval(f, x0+(h/2), y0+(k2/2),lam);     k4 = h*feval(f, x0+h, y0+k3,lam);       a=y0(1);       y0 = y0 + (k1 + 2*k2 + 2*k3 + k4)/6;     x0 = x0 + h;     values(i+1)=y0(1);       if sgn(a)~=sgn(y0(1))         spr = spr +1;     end       %hold on     %plot(x0,y0(1),'-x')     %plot(x0,y0(3),'-*')   end
I made also program for searching interval where n-th eigenvalue is:
Code:

function root = BISECTION(f,xa,y0,xb,lam,N,c,d,whichLambda)     root=(c+d)/2; while (root~=c) & (root~=d)       if whichLambda>SIGN(f,xa,root,lam,N)           c=root;     else         d=root;     end     root=(c+d)/2; end
SIGN counts number of sign changes in interval [xa, root]

And here is problem...newton iteration. My program:

Code:

function lam = fN(f, x0, y0, xb, lam, N, n) %lam = fN(f, x0, y0, xb, lam, N, n) %n = number of steps for Newton's iteration %N = number of steps for RK4   lam1=16.1;   for i = 1:n       [x1, y1, spr, vr] = RungeKutta4(f, x0, y0, xb, lam, N);       lam1 = lam1 - (y1(1))/(y1(3));   end
I'm not getting any good results for eigenvalue (lam1).

Thanks for any help...

Igor

Are you sure that you are supposed to be using nummerical integration on this problem?

.
• Jan 24th 2009, 08:31 AM
IgorMele
I have to do this in the second part. In the first part I used divided differences (approximatios for second and first derivations) and use eig(M), where M is tridiagonal matrix.

Ly = y'' + p(x)y' + q(x)y = lambda*y

My program:
Code:

function [lv,a,b,h2] = eigenvalues(p,q,a,b,n) %n = size of matrix h = (b-a)/n; x = a : h : b; h2 = h*h; qv = feval(q , x(2:end-1)); pv = (h/2)*feval(p , x(2:end-1)); Ma = diag(-2+h2*qv) + diag(1+pv(1:end-1),1) + diag(1-pv(2:end),-1); mia = eig(Ma); EigenValues = mia/h2;
Problem with this method is accuracy of higher eigenvalues. For example:
M = 1000×1000
100th eigen value supposed to be 10 000, matlab result: 9918.0234011
So I have to figure out how to calculate a specific eigenvalue with greater accuracy...
Here are my results:
http://www.shrani.si/t/1N/wE/17h0ThwZ/eig.jpg
• Jan 25th 2009, 10:03 AM
IgorMele
Hi!
I'm glad to report that my program is now working. It's basically the same code except:
Code:

lam = lam - (y1(1))/(y1(3));
The reason why I wasn't getting good results for $\lambda$ is too small number of steps in RK4 method.

Here are some results:
M = 1000 × 1000
33th eigen value supposed to be 33*33 = 1089,
matlab result: 1088.024968452657
Newton's method(3000 steps of RK4 from $[0,\pi]$):
1089.000025873885
Me very happy (Happy)