[Octave-bug-tracker] [bug #44004] eig and eigs routine

[Sebastian Glane]

Follow-up Comment #2, bug #44004 (project octave):

I modified the code to ensure convergence (flag -- fl) and to check if the
eigenvalues are real. If eigenvalues may be complex the sort command can cause
the incorrect results.

% solves the Eigenvalue Problem 
%   -y''= lambda^2 y on the interval (a,b)
%
% with mixed boundary conditions
%   -- dirichlet boundary condition y(a)=0
%   -- robin boundary condition   y'(b)=y(b)
clear all;
close all;
p=6;
%% number of points
n=2^p;
%% grid space
h=1/2^p;
%% numbers of eigenvalues to find
kmax=8;
%% assemble symmetric 2nd finite difference matrix
Lh=-gallery('tridiag',n,1,-2,1)./(h^2);
%% apply robin boundary in last line
Lh(n,:)=0.0;
Lh(n,n-1)=-2/h^2;
Lh(n,n)=-1/h^2*(2*h-2);
%% options for eigs
% opts.tol=1e-10;
opts.p=max(floor(n/2),2*kmax);
comp=0;
while comp<10 % do the same thing 10 times
        %% use eig to get all eigenvalues
        [~,LambdaEig]=eig(Lh);
        if isreal(diag(LambdaEig))
                LambdaEig=sort(diag(LambdaEig));
                lambdaEig=LambdaEig(1:kmax);
        else
                warning('eigenvalues of eig are not real')
        end
        %% use eigs routine to get 1st kmax eigenvalues of smallest magnitude
        [~,LambdaEigs,fl]=eigs(Lh,kmax,'sr',opts);
        if isreal(diag(LambdaEigs))
                LambdaEig=sort(diag(LambdaEig));
                lambdaEigs=sort(diag(LambdaEigs));
        else
                warning('eigenvalues of eigs are not real')
        end
        %% compute difference of 1st kmax eigenvalues of smallest magnitude
        errmax=max(abs(lambdaEig-lambdaEigs));
        disp(['eigs flag: ' num2str(fl) ', maximum error of eigenvalues: '...
                                num2str(errmax)])
        comp=comp+1;
end

    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?44004>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/