jacobi.m

来自「自己编写的程序」· M 代码 · 共 43 行

function jacobi(A,e,N)
k=1;
while k<=N
n=length(A);
b=abs(A-tril(A));
bb=max(max(b));
for i=1:n
    for j=i:n
        if b(i,j)==bb 
            t=i;
            l=j;
            break;
        end 
    end
end

ctg2=(A(t,t)-A(l,l))/(2*A(t,l));
s1=sin(acot(ctg2)/2);
c1=cos(acot(ctg2)/2);
V=zeros(n,n);
for p=1:n
    for q=1:n
        if p==q
            V(p,q)=1;
        else
            V(p,q)=0;
        end
    end
end
V(t,t)=c1;
V(l,l)=c1;
V(t,l)=-s1;
V(l,t)=s1;
A1=V'*A*V;
if abs(sum(sum(A1-diag(diag(A1)))))<=e
    A1
    break;
else
    k=k+1;
    A=A1;
end
end
%A=[2 -1 0;-1 2 -1;0 -1 2]