Thread: [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

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?

.

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:

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 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 ):
1089.000025873885
Me very happy