Results 1 to 4 of 4

Math Help - [MATLAB] Searching for specific eigenvalues with Newton iteration

  1. #1
    Newbie
    Joined
    Jan 2009
    Posts
    3

    Post [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
    Follow Math Help Forum on Facebook and Google+

  2. #2
    Member
    Joined
    May 2006
    Posts
    244
    Quote Originally Posted by IgorMele View Post
    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?

    .
    Follow Math Help Forum on Facebook and Google+

  3. #3
    Newbie
    Joined
    Jan 2009
    Posts
    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:
    Follow Math Help Forum on Facebook and Google+

  4. #4
    Newbie
    Joined
    Jan 2009
    Posts
    3

    Smile

    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
    Follow Math Help Forum on Facebook and Google+

Similar Math Help Forum Discussions

  1. [SOLVED] Numerical analysis, Newton's iteration
    Posted in the Differential Geometry Forum
    Replies: 8
    Last Post: March 27th 2011, 07:46 AM
  2. Replies: 1
    Last Post: February 7th 2010, 09:42 AM
  3. Newton-Raphson Iteration sequence
    Posted in the Calculus Forum
    Replies: 0
    Last Post: November 27th 2008, 12:48 PM
  4. Newton iteration; starting value?
    Posted in the Calculus Forum
    Replies: 1
    Last Post: October 19th 2008, 05:24 AM
  5. Iteration in MATLAB
    Posted in the Math Software Forum
    Replies: 4
    Last Post: June 20th 2008, 03:13 PM

Search Tags


/mathhelpforum @mathhelpforum