值分析第四版 颜庆津 计算实习题P238 结果：

%function [lam]=jingmi(ep) %反幂法求矩阵按模最小特征值 %2015.11.8 密密编写  (*^__^*) 嘻嘻……     function [lam]=jingmi(ep)%lam是A的按模最小的特征值 %初始化    n=501;    a=1:n;%存储A的主对角线元素    for i=1:n        a(i)=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i);    end    b=0.16;    c=-0.064;    u0=ones(n,1);%迭代初始向量    u1=ones(n,1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%第一次迭代            m=sqrt(u0'*u0);    y=u0/m; %以下代码功能：三角分解法，从线性方程组A*u1=y求解u1 b0=y;    L=zeros(3,n);    U=zeros(3,n);        %L中的主对角线元素均为1    for i=1:n         L(1,i)=1;    end %%%%%%%%%%%%%%%%%%%%%%doolittle分解,第一部分A=LU    for k=1:n               for j=k:min(k+2,n) %计算 U(k,j)            sum=0; %临时和            M=[1,k-2,j-2];%存储t的初始量备选值，t取M中的最大值            for t=max(M):k-1                 sum=sum+L(k-t+1,t)*U(t-j+3,j);            end           %%%%%%以下程序功能相当于语句U(k,j)=A(k,j)-sum; %计算L的第k列元素             if k==j                 U(k-j+3,j)=a(k)-sum;            else if j==k+1                     U(k-j+3,j)=b-sum;                 else if j==k+2                        U(k-j+3,j)=c-sum;                     end                 end           end           %%%%%%以上程序功能相当于语句U(k,j)=A(k,j)-sum; %计算L的第k列元素        end               for i=k+1:min(k+2,n)%计算 L(i,k)            sum=0; %临时和            M=[1,i-2,k-2];%存储t的初始量备选值，t取M中的最大值            for t=max(M):k-1                 sum=sum+L(i-t+1,t)*U(t-k+3,k);            end                      %%%%%%以下程序功能相当于语句  L(i,k)=(A(i,k)-sum)/U(k,k);%计算U的第k行元素            if k==i-1                 L(i-k+1,k)=(b-sum)/U(k-k+3,k);            else if k==i-2                    L(i-k+1,k)=(c-sum)/U(k-k+3,k);                 end            end            %%%%%%以上程序功能相当于语句L(i,k)=(A(i,k)-sum)/U(k,k);%计算U的第k行元素        end    end %%%%%%%%%%%%%%%%%%%%%%第二部分，求解Ly0=b1和Ux0=y0    y0=zeros(n,1);    x0=zeros(n,1);    y0(1)=b0(1);    for i=2:n %解y0        sum=0;        for t=max(1,i-2):i-1            sum=sum+L(i-t+1,t)*y0(t);        end        y0(i)=b0(i)-sum;    end    x0(n)=y0(n)/U(n-n+3,n);    for i=n-1:-1:1  %解x0        sum=0;        for t=i+1:min(i+2,n)            sum=sum+U(i-t+3,t)*x0(t);        end        x0(i)=(y0(i)-sum)/U(i-i+3,i);    end    %%%%%%%%%%%%%%%%%%%%doolittle三角分解求解线性方程组Ax=b，结束。    u1=x0;%返回结果 %以上代码功能：三角分解法，从线性方程组A*u1=y求解u1        x_prior=y'*u1;    err=1;    u0=u1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%循环迭代开始  while(err>ep)    m=sqrt(u0'*u0);    y=u0/m; %以下代码功能：三角分解法，从线性方程组A*u1=y求解u1  b0=y;    L=zeros(3,n);    U=zeros(3,n);        %L中的主对角线元素均为1    for i=1:n        L(i-i+1,i)=1;    end %%%%%%%%%%%%%%%%%%%%%%doolittle分解,第一部分A=LU    for k=1:n               for j=k:min(k+2,n) %计算 U(k,j)            sum=0; %临时和            M=[1,k-2,j-2];%存储t的初始量备选值，t取M中的最大值            for t=max(M):k-1                 sum=sum+L(k-t+1,t)*U(t-j+3,j);            end           %%%%%%以下程序功能相当于语句U(k,j)=A(k,j)-sum; %计算L的第k列元素             if k==j                 U(k-j+3,j)=a(k)-sum;            else if j==k+1                     U(k-j+3,j)=b-sum;                 else if j==k+2                         U(k-j+3,j)=c-sum;                     end                 end           end           %%%%%%以上程序功能相当于语句U(k,j)=A(k,j)-sum; %计算L的第k列元素        end               for i=k+1:min(k+2,n)%计算 L(i,k)            sum=0; %临时和            M=[1,i-2,k-2];%存储t的初始量备选值，t取M中的最大值            for t=max(M):k-1                 sum=sum+L(i-t+1,t)*U(t-k+3,k);            end                      %%%%%%以下程序功能相当于语句  L(i,k)=(A(i,k)-sum)/U(k,k);%计算U的第k行元素            if k==i-1                 L(i-k+1,k)=(b-sum)/U(k-k+3,k);            else if k==i-2                    L(i-k+1,k)=(c-sum)/U(k-k+3,k);                 end            end            %%%%%%以上程序功能相当于语句L(i,k)=(A(i,k)-sum)/U(k,k);%计算U的第k行元素        end    end %%%%%%%%%%%%%%%%%%%%%%第二部分，求解Ly0=b1和Ux0=y0    y0=zeros(n,1);    x0=zeros(n,1);    y0(1)=b0(1);    for i=2:n %解y0        sum=0;        for t=max(1,i-2):i-1            sum=sum+L(i-t+1,t)*y0(t);        end        y0(i)=b0(i)-sum;    end    x0(n)=y0(n)/U(n-n+3,n);    for i=n-1:-1:1  %解x0        sum=0;        for t=i+1:min(i+2,n)            sum=sum+U(i-t+3,t)*x0(t);        end        x0(i)=(y0(i)-sum)/U(i-i+3,i);    end    %%%%%%%%%%%%%%%%%%%%doolittle三角分解求解线性方程组Ax=b，结束。    u1=x0;%返回结果 %以上代码功能：三角分解法，从线性方程组A*u1=y求解u1    lam=y'*u1;    err=abs(lam-x_prior)/abs(lam);    u0=u1;    x_prior=lam;   end    lam=y'*u1;    err=abs(lam-x_prior)/abs(lam);    u0=u1;    x_prior=lam; end lam=1/lam; end