C..4 シフト法

% eigen3.m --- shift method
format long
a=[1 1 1;1 2 2;1 2 3]
eigen=eig(a)
x=ones(3,1);
maxiter=10;
r=zeros(maxiter,1);
shift=0.6;
% we cannot use the factorization of a, can we?
[L u p]=lu(a-shift*eye(3,3));
for k=1:maxiter
    y=u\(L\(p*x));
    x=y/norm(y);
    r(k)=x'*(a*x);
end
r
eigen=sort(eigen);
error=r-eigen(2)
clf
semilogy(1:maxiter,abs(error))