1、第一題 關于舍入誤差累計的效果模擬。x*是在（0，1）上服從均勻分布的隨機數，對x*取5位有效數字得到x,將產生的k（比如取為10000）個x*（x）相加得到X*（X），研究X*-X的分布情況以及X*-X和k的關系。將得到的結果用圖形表示出來。k=1;kk=ones(1,1000);xc=ones(1,1000);while k<=1000r=rand(1,k);v=vpa(r,5);x0=sum(r);x1=sum(v);cha=x0-x1;xc(k)=cha;kk(k)=k;k=k+1;end;plot(kk,xc,'rh') 可以看出隨著處理的數據的增大誤

2、差也越來越大，但是還分布在橫軸的兩側第二題 研究產生各種特定矩陣的方法（階數在10-100），比如對稱陣，三對角陣，正定矩陣，正交陣，對角占優矩陣，說明如何生成。（1）生成對稱矩陣k=10;a=rand(10,10);for i=1:k %先生成一個下三角矩陣; for j=1:k if j>i a(i,j)=0; end; end;end;aa=a' %加上他本身的轉置;g=a+aag = 1.9760 0.0377 0.8852 0.9133 0.7962 0.0987 0.2619 0.3354 0.6797 0.1366 0.0377 0.2135 0.6538 0.49

3、42 0.7791 0.7150 0.9037 0.8909 0.3342 0.6987 0.8852 0.6538 1.4881 0.5000 0.4799 0.9047 0.6099 0.6177 0.8594 0.8055 0.9133 0.4942 0.5000 1.7730 0.0287 0.4899 0.1679 0.9787 0.7127 0.5005 0.7962 0.7791 0.4799 0.0287 0.1429 0.5216 0.0967 0.8181 0.8175 0.7224 0.0987 0.7150 0.9047 0.4899 0.5216 1.6007 0.4

4、538 0.4324 0.8253 0.0835 0.2619 0.9037 0.6099 0.1679 0.0967 0.4538 0.7985 0.5269 0.4168 0.6569 0.3354 0.8909 0.6177 0.9787 0.8181 0.4324 0.5269 0.7448 0.1981 0.4897 0.6797 0.3342 0.8594 0.7127 0.8175 0.8253 0.4168 0.1981 1.8855 0.4177 0.1366 0.6987 0.8055 0.5005 0.7224 0.0835 0.6569 0.4897 0.4177 1.

5、9982（2）三對角矩陣k=zeros(5,5);a=rand(1,5);b=rand(1,4);c=rand(1,4);a=diag(a);b=diag(b);c=diag(c);k(1:4,2:5)=c;aa=k; %將得到的結果賦值給其他變量 因為此時的K已經不是以前的了必須再次初始化 才能不會影響以后的操作（如下一步）k=zeros(5,5);k(2:5,1:4)=b;bb=k;aa+bb+a %將上述的三個矩陣復合得到三對角;ans = 0.5005 0.8175 0 0 0 0.0714 0.4711 0.7224 0 0 0 0.5216 0.0596 0.1499 0 0 0

6、0.0967 0.6820 0.6596 0 0 0 0.8181 0.0424也是三對角k=6;a=zeros(6,6);a(1,1:2)=rand(1,2);a(6,5:6)=rand(1,2);i=2; %i不能忘了初始化while i>1&i<6a(i,(i-1):(i+1)=rand(1,3);i=i+1;end;aa = 0.1389 0.6963 0 0 0 0 0.5303 0.8611 0.4849 0 0 0 0 0.3935 0.6714 0.7413 0 0 0 0 0.5201 0.3477 0.1500 0 0 0 0 0.5861 0.2621

7、 0.0445 0 0 0 0 0.0938 0.5254（3）正定矩陣a=rand(5,5);b=a'c=a*b c = 0.3505 0.6554 0.6104 0.5783 0.7651 0.6554 1.6287 1.4020 0.9147 1.5272 0.6104 1.4020 1.4758 1.1823 1.7376 0.5783 0.9147 1.1823 1.4337 1.7436 0.7651 1.5272 1.7376 1.7436 2.3737矩陣的轉置和他本身的乘積是正定的（4）正交矩陣 隨機生成一個方陣，然后利用schmidt正交化對列向量進行正交化單位化得

8、到一個正交矩陣；a=rand(5,5);a=vpa(a,7);b=zeros(5,5);b(:,1)=a(:,1);for i=2:5sum=zeros(5,1); for j=1:(i-1) sum=sum+(dot(a(:,i),b(:,j)/dot(b(:,j),b(:,j)*b(:,j); end; b(:,i)=a(:,i)-sum;end; %完成對矩陣的正交化schmidt;for k=1:5b(:,k)=b(:,k)/sqrt(dot(b(:,k),b(:,k); end; %完成對矩陣的單位化;b=vpa(b,7)b*b' b = 0.4662487, -0.5697

9、477, 0.4464938, 0.5031049, 0.07435341 0.4330088, 0.5708757, -0.2023893, 0.4915742, -0.451661 0.6966242, -0.2311083, -0.434368, -0.5221337, 0.002112758 0.2120276, 0.3050953, 0.7502452, -0.4773735, -0.266848 0.2547048, 0.4505489, 0.09021349, 0.06878292, 0.8480929 ans = 1.000000000000000658769429433138

10、, 0.00000000000000020501909475855146103181737845679, -0.00000000000000021440423350092585100221961008361, 0.0000000000000005427542082813480164677324711907, -0.00000000000000011700162166438432633834107099474 0.00000000000000020501909475855146103181737845679, 1.0000000000000003207231381623452, 0.000000

11、00000000058082238049994760442870183594865, -0.00000000000000073536380504684284367618389622692, 0.00000000000000040560052001944880242963606302516 -0.00000000000000021440423350092585100221961008361, 0.00000000000000058082238049994760442870183594865, 1.0000000000000016447931562774053, 0.000000000000000

12、63268757033543168173701894986961, -0.00000000000000013226852520735450182808565768223 0.0000000000000005427542082813480164677324711907, -0.00000000000000073536380504684284367618389622692, 0.00000000000000063268757033543168173701894986961, 0.99999999999999936362009174092183, 0.000000000000000255697544

13、97442931226184064248466 -0.00000000000000011700162166438432633834107099474, 0.00000000000000040560052001944880242963606302516, -0.00000000000000013226852520735450182808565768223, 0.00000000000000025569754497442931226184064248466, 0.99999999999999902311091880723728（5）對角占優矩陣隨機生成一個矩陣，將每行的元素求和然后復賦值給對角元素

14、，得到對角占優矩陣；a=4*rand(10,10);for i=1:6a(i,i)=sum(abs(a(i,1:6); %對角占優矩陣式對角線元素的絕對值之和大于其他元素的絕對值之和，所以這里取了絕對值;endaa = 12.3828 1.3800 0.2741 3.0666 3.6312 1.7501 0.9074 2.7388 0.3920 2.7565 1.3859 8.8416 0.8722 0.0873 3.8983 1.2834 1.7841 1.7724 0.9827 2.8715 2.2300 3.7099 12.0057 1.5724 0.4795 0.5362 1.0649

15、 1.7427 2.4629 2.2361 1.1991 3.0243 1.6569 9.5047 2.0759 0.5383 1.8364 3.1721 1.2198 2.1334 0.6363 1.1530 2.6448 0.8169 11.7626 3.2238 1.7316 3.2622 3.0679 3.5029 2.6610 2.4247 3.1339 2.6491 2.5481 15.5160 1.0385 3.0085 1.0689 1.5724 2.7368 3.0641 0.9915 3.6590 3.8156 3.7770 0.5349 3.1570 0.1580 1.8

16、323 3.1696 3.3846 2.2177 0.0276 3.7877 3.9534 1.6769 2.0051 1.1862 0.8330 1.3945 3.6079 0.9183 2.9857 3.8664 1.6396 2.0274 2.2207 2.2255 3.02911.0003 2.3828 0.0277 3.1987 0.2694 1.4848 1.2973 2.5230 3.8763 2.1868下面也是三對角占優k=6;a=zeros(6,6);a(1,1:2)=rand(1,2);a(6,5:6)=rand(1,2);i=2; %i不能忘了初始化哦親while i&

17、gt;1&i<6a(i,(i-1):(i+1)=rand(1,3);i=i+1;end;for i=1:6 a(i,i)=sum(a(i,:);end; %在三對角矩陣的基礎上將每一行的數相加然后復制給對角元素 則得到的矩陣式對角占優的;aa = 1.0292 0.8620 0 0 0 0 0.8843 1.6271 0.1548 0 0 0 0 0.1999 1.3555 0.7487 0 0 0 0 0.8256 1.9341 0.3185 0 0 0 0 0.5341 0.7357 0.1117 0 0 0 0 0.9899 1.5043第三題 比較高斯消去法和帶有列主元的

18、高斯消去法的效果。要求構造一些方程組，方程組的精確解已知，用上面兩種方法分別求解，說明選主元的效果。高斯消去法 A=1 2 3 1;2 7 5 6;1 4 9 -3; %增廣矩陣;x=zeros(1,3);m,n=size(A);for i=1:(m-1) for j=(i+1):m A(j,:)=A(j,:)-A(i,:)*A(j,i)/A(i,i); %消去;end endA;x(m)=A(m,n)/A(m,m); %回代 未知數的最后一個分量單獨計算出其他的如下;for i=(m-1):-1:1 x(i)=(A(i,n)-A(i,i+1:m)*x(i+1:m)')/A(i,i);

19、 end xx = 2 1 -1 高斯列主元法A=1 2 3 1;2 7 5 6;1 4 9 -3; %增廣矩陣;m,n=size(A);for i=1:(m-1) temp=max(abs(A(i:m,i); %當前處理的矩陣的第一列的絕對值最大的元素; a,b=find(abs(A(i:m,i)=temp); %找到最大元素所在的位置; tempo=A(a(1)+i-1,:); A(a(1)+i-1,:)=A(i,:); A(i,:)=tempo; %交換兩行;for j=(i+1):m A(j,:)=A(j,:)-A(i,:)*A(j,i)/A(i,i); %消去; end A;endx

20、(m)=A(m,n)/A(m,m); %回代;for i=(m-1):-1:1 x(i)=(A(i,n)-A(i,i+1:m)*x(i+1:m)')/A(i,i); end xA = 2 7 5 6 0 -3/2 1/2 -2 0 1/2 13/2 -6 A = 2 7 5 6 0 -3/2 1/2 -2 0 0 20/3 -20/3 x = 2 1 -1 下面構造一些已知解的矩陣，然后分別用以上兩種方法進行求解，比較它們運行的結果：高斯消去A=rand(500,500);yizhij=rand(500,1);b=A*yizhij;A=A,b; %增廣矩陣;x=zeros(1,500)

21、;m,n=size(A);for i=1:(m-1) for j=(i+1):m A(j,:)=A(j,:)-A(i,:)*A(j,i)/A(i,i); %消去;end endA;x(m)=A(m,n)/A(m,m); %回代 未知數的最后一個分量單獨計算出其他的如下;for i=(m-1):-1:1 x(i)=(A(i,n)-A(i,i+1:m)*x(i+1:m)')/A(i,i); end dot(x'-yizhij,x'-yizhij)ans = 4.8354e-020高斯列主元法A=rand(500,500);yizhij=rand(500,1);b=A*yiz

22、hij;A=A,b; %增廣矩陣; m,n=size(A);for i=1:(m-1) temp=max(abs(A(i:m,i); %當前處理的矩陣的第一列的絕對值最大的元素; a,b=find(abs(A(i:m,i)=temp); %找到最大元素所在的位置; tempo=A(a(1)+i-1,:); A(a(1)+i-1,:)=A(i,:); A(i,:)=tempo; %交換兩行;for j=(i+1):m A(j,:)=A(j,:)-A(i,:)*A(j,i)/A(i,i); %消去; end A;endx(m)=A(m,n)/A(m,m); %回代;for i=(m-1):-1:1

23、 x(i)=(A(i,n)-A(i,i+1:m)*x(i+1:m)')/A(i,i); end dot(x'-yizhij,x'-yizhij)ans = 6.1236e-023可以看出運用高斯列主元法求解比高斯消去法更加精確。第四題采用緊湊格式實現矩陣的LU分解，并使用這個方法求解線性方程組。A=2 4 2 6;4 9 6 15;2 6 9 18;6 15 18 40;N = size(A);n = N(1);L = zeros(4,4);U = eye(4,4); %U的對角元素為1;L(1:4,1) = A(1:4,1); %L的第一列;U(1,1:4) = A(

24、1,1:4)/L(1,1); %U的第一行;for k=2:4 for i=k:4 L(i,k) = A(i,k)-L(i,1:(k-1)*U(1:(k-1),k); %L的第k列; end for j=(k+1):4 U(k,j) = (A(k,j)-L(k,1:(k-1)*U(1:(k-1),j)/L(k,k); %U的第k行; endendLUb=9 23 22 47'y = inv(L)*b;x = inv(U)*y %求解方程;L = 2 0 0 0 4 1 0 0 2 2 3 0 6 3 6 1 U = 1 2 1 3 0 1 2 3 0 0 1 2 0 0 0 1 x =

25、 1/2 2 3 -1 第五題 設計對稱正定線性方程組，使用平方根法計算方程組的解。A=2 4 2 6;4 9 6 15;2 6 9 18;6 15 18 40;b=9 23 22 47'n=length(b);%方程個數nL=zeros(n,n);L(1,1)=sqrt(A(1,1);L(2:n,1)=A(2:n,1)/L(1,1);for j=2:n-1 L(j,j)=sqrt(A(j,j)-sum(L(j,1:j-1).2); for i=j+1:n L(i,j)=(A(i,j)-sum(L(i,1:j-1).*L(j,1:j-1)/L(j,j); endendL(n,n)=sq

26、rt(A(n,n)-sum(L(n,1:n-1).2);L y = inv(L)*b;x = inv(L')*y %求解方程;L = 1393/985 0 0 0 3363/1189 1 0 0 1393/985 2 1351/780 0 4756/1121 3 1351/390 1 x = 1/2 2 3 -1 第六題 追趕法求解三對角線性方程組三對角線性方程組的追趕法編程，使得算法計算次數達到最小。a=2.0 -1 0 0 0 ;-1 2 -1 0 0;0 -1 2 -1 0;0 0 -1 2 -1;0 0 0 -1 2;b=1 0 0 0 0;y=zeros(1,5);x=zer

27、os(1,5);L=zeros(5,5);U=zeros(5,5);U(1,2)=a(1,2)/a(1,1);for i=2:4 %利用書上給出的步驟編寫; U(i,(i+1)=a(i,(i+1)/(a(i,i)-a(i,(i-1)*U(i-1),i);end;y(1)=b(1)/a(1,1);for i=2:5 y(i)=(b(i)-a(i,(i-1)*y(i-1)/(a(i,i)-a(i,(i-1)*U(i-1),i);end;x(5)=y(5);i=4;while i>=1 x(i)=y(i)-U(i,(i+1)*x(i+1); i=i-1;end;yx=vpa(x,4)y = 1

28、/2 1/3 1/4 1/5 1/6 x = 0.8333, 0.6667, 0.5, 0.3333, 0.1667第七題 調用Matlab的cond()函數分析隨機生成的非奇異陣的的階數和條件數之間的關系i=11;kk=zeros(1,90);con=zeros(1,90);while i<=100 kk(i-10)=i; k=1; chengji=1; while k<=50 a=rand(i,i); %生成i階的方陣，對么一個方陣作50次的求解條件數的運算，最后求取平均值; c=cond(a); chengji=chengji*c; %把每次一運行出來的結果乘在一起，最后取幾

29、何平均數; k=k+1; end; c=chengji.(1/50); con(i-10)=c; i=i+1;end;plot(kk,con,'rh')hold on;p=polyfit(kk,con,2);con=polyval(p,kk);plot(kk,con)可以看出隨著非奇異矩陣階數的不斷增大條件數呈現多項式函數增大第八題 編寫線性方程組的Jacobi方法和高斯賽德爾迭代法的程序。使用圖形反映譜半徑和殘差之間的關系。Jacobi方法高斯賽德爾迭代法A=8 -3 2;4 11 -1;6 3 12;b=20 33 36'e=0.000001;n=max(size(

30、A);for i=1:n if A(i,i)=0 '對角元素為0不能求解' return; end;end;x=ones(n,1);k=0; %迭代次數的初始化;mk=50; %迭代次數的最大值;r=1;%前后項之差的無窮范數;while k<=mk&r>e x0=x; for i=1:n sum=0; for j=1:i-1 sum=sum+A(i,j)*x0(j); end; for j=i+1:n sum=sum+A(i,j)*x0(j); end; x(i)=(b(i)-sum)/A(i,i); end; %以上是迭代公式的等價翻譯;r=norm(x

31、-x0,inf);k=k+1;end;if k>mk '迭代失敗'else'迭代成功'end;xkans =迭代成功x = 3.0000 2.0000 1.0000k = 15A=8 -3 2;4 11 -1;6 3 12;b=20 33 36'e=0.000001;n=max(size(A);for i=1:n if A(i,i)=0 '對角元素為0不能求解' return; end;end;x=ones(n,1);k=0; %迭代次數的初始化;mk=50; %迭代次數的最大值;r=1;%前后項之差的無窮范數;while k<

32、;=mk&r>e x0=x; for i=1:n sum=0; for j=1:i-1 sum=sum+A(i,j)*x(j); end; for j=i+1:n sum=sum+A(i,j)*x0(j); end; x(i)=(b(i)-sum)/A(i,i); end; %以上是迭代公式的等價翻譯;r=norm(x-x0,inf);k=k+1;end;if k>mk '迭代失敗'else'迭代成功'end;xkans =迭代成功x = 3.0000 2.0000 1.0000k = 8其中高斯賽德爾迭代法只是在Jacobi方法在一個地方的

33、更改。并且可以看出高斯賽德爾迭代法比Jacobi方法收斂更快。第九題 編寫線性方程組的SOR方法的程序，對一個給定的矩陣，選擇不同的不同的松弛因子，研究松弛因子和迭代矩陣的譜半徑之間的關系。w=1w=1.5A=8 -3 2;4 11 -1;6 3 12;b=20 33 36'e=0.000001;w=1;n=max(size(A);for i=1:n if A(i,i)=0 '對角元素為0不能求解' return; end;end;x=ones(n,1);k=0; %迭代次數的初始化;mk=100; %迭代次數的最大值;r=1;%前后項之差的無窮范數;while k&l

34、t;=mk&r>e x0=x; for i=1:n sum=0; for j=1:i-1 sum=sum+A(i,j)*x(j); end; for j=i+1:n sum=sum+A(i,j)*x0(j); end; x(i)=(1-w)*x0(i)+w/A(i,i)*(b(i)-sum); end; %以上是迭代公式的等價翻譯;r=norm(x-x0,inf);k=k+1;end;if k>mk '迭代失敗'else'迭代成功'end;wxkans =迭代成功w = 1x = 3.0000 2.0000 1.0000k = 8A=8 -3

35、 2;4 11 -1;6 3 12;b=20 33 36'e=0.000001;w=1.5;n=max(size(A);for i=1:n if A(i,i)=0 '對角元素為0不能求解' return; end;end;x=ones(n,1);k=0; %迭代次數的初始化;mk=100; %迭代次數的最大值;r=1;%前后項之差的無窮范數;while k<=mk&r>e x0=x; for i=1:n sum=0; for j=1:i-1 sum=sum+A(i,j)*x(j); end; for j=i+1:n sum=sum+A(i,j)*x0

36、(j); end; x(i)=(1-w)*x0(i)+w/A(i,i)*(b(i)-sum); end; %以上是迭代公式的等價翻譯;r=norm(x-x0,inf);k=k+1;end;if k>mk '迭代失敗'else'迭代成功'end;wxkans =迭代成功w = 1.5000x = 3.0000 2.0000 1.0000k = 51松弛因子和迭代矩陣的譜半徑之間的關系A=8 -3 2;4 11 -1;6 3 12;for i=1:3 for j=1:3 if j>=i A(i,j)=0; end; end;end;L=A;A=8 -3

37、2;4 11 -1;6 3 12;for i=1:3 for j=1:3 if j<=i A(i,j)=0; end; end;end;U=A;A=8 -3 2;4 11 -1;6 3 12;D=A-L-U;L,D,U,pbj=zeros(1,199);kk=zeros(1,199);w=0.01;k=1; %這里用K是為了防止下面出現錯誤;while w<2 B=inv(D+w*L)*(1-w)*D-w*U); p=vrho(B); pbj(k)=p; %記錄下譜半徑; kk(k)=w; %記錄下松弛因子; k=k+1; w=w+0.01;end;plot(kk,pbj,'

38、;rh') %畫圖; hold on;p=polyfit(kk,pbj,2)pbj=polyval(p,kk);plot(kk,pbj,'b') %畫圖;為什么1.5的哪個位置是一個轉折點呢第十題 編寫線性方程組的最速下降法的程序，并研究對于給定的矩陣A（給定特征值），設定矩陣A的所有特征值，研究殘差的范數和迭代次數k之間的關系，說明最速下降法的缺點在哪里。同時研究特征值的變化是否影響迭代收斂的速度。最速下降法求解線性方程組的程序和運行結果如下:e=0.0000001; %給定誤差限;A=4 -1 2;-1 5 3;2 3 6;b=12;10;18;x=1;1;1;k=

39、1; %初始化，是用來統計迭代次數的;r=b-A*x; %初始化r，并用于第一次的判斷“while norm(r)>e”;while norm(r)>er=b-A*x;arf=dot(r,r)/dot(A*r,r); %步長的計算公式;y=x+arf*r; %得到下一個迭代點;r=b-A*y; %得到新的剩余向量;x=y; k=k+1; %得帶次數加1;end;k %輸出迭代次數和最后一次迭代得到的x;xk = 76x = 3.0000 2.0000 1.0000研究殘差的范數和迭代次數k之間的關系(二次擬合的效果)e=0.0000001; %給定誤差限;A=4 -1 2;-1 5

40、 3;2 3 6;b=12;10;18;x=1;1;1;kk=zeros(1,1000);rr=zeros(1,1000);k=1; %初始化，是用來統計迭代次數的;r=b-A*x; %初始化r，并用于第一次的判斷“while norm(r)>e”;while norm(r)>er=b-A*x;arf=dot(r,r)/dot(A*r,r); %步長的計算公式;y=x+arf*r; %得到下一個迭代點;r=b-A*y; %得到新的剩余向量;rr(k)=norm(r);x=y; kk(k)=k;k=k+1; %得帶次數加1;end;kkk=kk(1,1:k);rr=rr(1,1:k)

41、;p=polyfit(kk,rr,2)rr=polyval(p,kk);plot(kk,rr)k = 76p = 0.0005 -0.0443 0.9259研究殘差的范數和迭代次數k之間的關系(三次擬合的效果)研究殘差的范數和迭代次數k之間的關系(四次擬合的效果)可以看出最速下降法收斂比較慢，事實上最速下降法只是線性收斂。當矩陣的特征值改變后：迭代次數變為kk=48;說明特征值的變化對收斂速度是有影響的，但是當特征值變化不大時對收斂速度的影響是比較小的第十一題 使用乘冪法計算方陣的主特征值和對應的特征向量。特別的要研究矩陣的特征值對于收斂速度的影響，分下面的情況：（1）主特征值是單根的情況，

42、（2）主特征值是重根的情況，（3）主特征值是共軛根的情況。（1）乘冪法計算方陣的主特征值和對應的特征向量A=1 -1 2;-2 0 5;6 -3 6;x0=1;1;1;eps = 1.0e-5;v = x0; M = 5000; m = 0; l = 0; %這是L;for(k=1:M) y = A*v; m = max(y); v = y/m; if(abs(m - l)<eps) %臨近的兩次得到的特征值的差的絕對值小于誤差限; l = m; s = k; else if(k=M) disp('迭代步數太多，收斂速度太慢！'); l = m; s = M; end e

43、ndendlvs迭代步數太多，收斂速度太慢！l = 5.0000v = 0.2778 0.8889 1.0000s = 5000矩陣的特征值對于收斂速度的影響，分下面的情況：（1）主特征值是單根的情況第十二題 對稱正定矩陣的加速迭代方法的驗證。第十三題 反冪法求矩陣的最小特征值的程序，為簡單起見，選擇矩陣是對稱正定的。在計算過程中，一定要使用LU分解的技術。第十四題 使用Householder方法將對稱矩陣化成三對角陣。A=6 2 3 1;2 5 4 8;3 4 9 1;1 8 1 7;B=zeros(m,n);m,n=size(A);i=m-2;k=1;for i=1:im,n=size(A);v1=A(1,2:n)' %依照筆記上的操作，編寫代碼如下;sgm=sqrt(dot(v1,v1);E=eye(n-1,1);v=(v1+sgm*E)./sqrt(dot(v1+sgm*E,v1+sgm*E);p1=eye(n-1)-2*(v*v');p=eye(n);p(2:m,2:n)=p1;A=p'*A*p;B(k:k+1,k:k+1)=A(k:k+1,k:k+1); %這里是為了將每次迭代產生的“A”的前兩行和前兩列