(1.11)-2023數(shù)值實驗-練習11

一、冪法（乘冪法）與反冪法第八章：特征值問題的計算方法1、冪法求按模最大的特征值和特征向量2、冪法的原點平移：加速收斂3、反冪法：求按模最小的特征值4、帶位移反冪法：求特征向量（高精度）二、Jacobi迭代法求對稱矩陣的全部特征值和特征向量5、Givens（旋轉）變換6、用Givens變換把向量的某些分量化為零7、用Givens變換把矩陣的某些元素化為零clear,clcA=rand(4);A=(A+A’)/2;n=length(A);k=0;whilenorm(tril(A,-1),'fro')>10^(-12)&k<10^4k=k+1;forp=1:n-1forq=p+1:nifA(p,q)~=0

ta=(A(q,q)-A(p,p))/(2*A(p,q));t=sign(ta)/(abs(ta)+sqrt(1+ta^2));c=1/sqrt(1+t^2);s=t*c;endfori=1:nif(i~=p&i~=q)u=c*A(i,p)-s*A(i,q);v=s*A(i,p)+c*A(i,q);A(i,p)=u;A(i,q)=v;A(p,i)=A(i,p);A(q,i)=A(i,q);endend

u=c^2*A(p,p)-2*s*c*A(p,q)+s^2*A(q,q);v=s^2*A(p,p)+2*s*c*A(p,q)+c^2*A(q,q);A(p,p)=u;A(q,q)=v;A(q,p)=0;A(p,q)=0;endendendA,k8、循環(huán)Jacobi方法求對稱矩陣的全部特征值和特征向量若：則先計算：再計算：最后形成：三、QR方法9、

Househlder變換關鍵步：令：如果再計算：則：其余分量

不唯一，也可以替換；為了唯一，此時，；省1個單元9、

Househlder變換function[u,beta]=house(x)n=length(x);x=x/norm(x,inf);de=x(2:n)'*x(2:n);u(1,1)=1;u(2:n,1)=x(2:n,1);ifde==0beta=0;elseaa=sqrt(x(1)^2+de);ifx(1)<=0u(1)=x(1)-aa;elseu(1)=-de/(x(1)+aa);endbeta=2*u(1)^2/(de+u(1)^2);u=u/u(1);end9、

Househlder變換若已知按模最大特征值及其特征向量：10、收縮技巧（降階方法）：求第二大特征值求一Householder矩陣H使得：則：即：則：7A=[210;131;014];lambda=4.7321;n=length(A);x=[0.2679,0.7321,1.0000]’;[u,beta]=house(x)H=eye(3)-beta*u*u’;A1=H*A*HB=A1(2:n,2:n);k=0;u0=ones(n-1,1);u=B*u0;[t,i]=max(abs(u));v=u(i);u=u/v;whilenorm(u-u0)>10^(-7)&k<10000u0=u;k=k+1;u=B*u0;[t,i]=max(abs(u));v=u(i);u=u/v;endu,v,k10、收縮技巧（降階方法）：求第二大特征值8A=[210;131;014];lambda=4.73205;n=length(A);x=[0.26795,0.73205,1.0000]’;B=A1(2:n,2:n);k=0;u0=ones(n-1,1);u=B*u0;[t,i]=max(abs(u));v=u(i);u=u/v;whilenorm(u-u0)>10^(-7)&k<10000u0=u;k=k+1;u=B*u0;[t,i]=max(abs(u));v=u(i);u=u/v;endu,v,k[u,beta]=house(x)%H=eye(3)-beta*u*u';B=H*A*Hw1=beta*(A*u);w2=w1-0.5*beta*(u’*w1)*u;A1=A-u*w2’-w2*u’;10*、收縮技巧（降階方法）：求第二大特征值clear,clcA=[1,3;3,5];k=0;whilenorm(tril(A,-1),’fro’)>10^(-6)&k<1000k=k+1;[Q,R]=qr(A);A=R*Q;endA,k11、基本的QR迭代clear,clcA=rand(5);%A=A’*A;n=length(A);fork=1:n-2[v,beta]=house(A(k+1:n,k));A(k+1:n,k:n)=(eye(n-k)-beta*v*v')*A(k+1:n,k:n);A(1:n,k+1:n)=A(1:n,k+1:n)*(eye(n-k)-beta*v*v');endA12、計算矩陣的上Hessenberg分解clear,clcA=rand(5);%A=A’*A;n=length(A);fork=1:n-2[v,beta]=house(A(k+1:n,k));A(k+1:n,k:n)=A(k+1:n,k:n)-beta*v*(v'*A(k+1:n,k:n));A(1:n,k+1:n)=A(1:n,k+1:n)-beta*(A(1:n,k+1:n)*v)*v';endA12*、計算矩陣的上Hessenberg分解一、Euler方法第九章：常微分方程初值問題數(shù)值解法1、向前Euler

公式2、向后Euler

公式3、梯形公式精確解：a=0;b=1;x(1)=a;y(1)=1;h=0.1;N=(b-a)/h;forn=1:Nx(n+1)=x(n)+h;y(n+1)=h*x(n)+(1-h)*y(n)+h;endplot(x,y,'ob');holdonforn=1:N

y(n+1)=(h*x(n+1)+y(n)+h)/(1+h);endplot(x,y,'*r');xx=a:0.01:b;yy=xx+exp(-xx);plot(xx,yy,‘-k’)1、向前和向后Euler法精確解：clear,clca=0;b=1;x(1)=a;y(1)=1;h=0.1;N=(b-a)/h;forn=1:Nx(n+1)=x(n)+h;y(n+1)=((1-h/2)