Thread: [MATLAB] Searching for specific eigenvalues with Newton iteration

1. [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:

then I found solution with RungeKutta method:

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

2. 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:

then I found solution with RungeKutta method:

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?

.

3. 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:

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