《MATLAB解线性方程组的直接方法.docx》由会员分享，可在线阅读，更多相关《MATLAB解线性方程组的直接方法.docx（13页珍藏版）》请在三一文库上搜索。
1、. 解线性方程组的直接方法在这章中我们要学习线性方程组的直接法，特别是适合用数学软件在计算机上求解的方法.3.1 方程组的逆矩阵解法及其MATLAB程序3.1.3 线性方程组有解的判定条件及其MATLAB程序判定线性方程组是否有解的MATLAB程序function RA,RB,n=jiepb(A,b)B=A b;n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA;if zhica>0,disp('请注意：因为RA=RB，所以此方程组无解.')returnendif RA=RB if RA=ndisp('请注意：因为RA=
2、RB=n，所以此方程组有唯一解.') else disp('请注意：因为RA=RB<n，所以此方程组有无穷多解.')endend例3.1.4 判断下列线性方程组解的情况.如果有唯一解，则用表 3-2方法求解.(1) (2) (3) (4) 解 在MATLAB工作窗口输入程序>> A=2 3 -1 5;3 1 2 -7;4 1 -3 6;1 -2 4 -7; b= 0; 0; 0; 0; RA,RB,n=jiepb(A,b)运行后输出结果为请注意：因为RA=RB=n，所以此方程组有唯一解.RA = 4,RB =4,n =4在MATLAB工作窗口输入>
3、;>X=Ab, 运行后输出结果为 X =(0 0 0 0).(2) 在MATLAB工作窗口输入程序>> A=3 4 -5 7;2 -3 3 -2;4 11 -13 16;7 -2 1 3;b= 0; 0; 0; 0;RA,RB,n=jiepb(A,b)运行后输出结果请注意：因为RA=RB<n，所以此方程组有无穷多解.RA =2,RB =2,n =4(3)在MATLAB工作窗口输入程序>> A=4 2 -1;3 -1 2;11 3 0; b=2;10;8; RA,RB,n=jiepb(A,B)运行后输出结果请注意：因为RA=RB，所以此方程组无解.RA =2,
4、RB =3,n =3(4)在MATLAB工作窗口输入程序>> A=2 1 -1 1;4 2 -2 1;2 1 -1 -1; b=1; 2; 1; RA,RB,n=jiepb(A,b)运行后输出结果请注意：因为RA=RB<n，所以此方程组有无穷多解.RA =2,RB =2,n =33.2 三角形方程组的解法及其MATLAB程序3.2.2 解三角形方程组的MATLAB程序解上三角形线性方程组的MATLAB程序function RA,RB,n,X=shangsan(A,b)B=A b; n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA;
5、if zhica>0,disp('请注意：因为RA=RB，所以此方程组无解.')returnendif RA=RB if RA=ndisp('请注意：因为RA=RB=n，所以此方程组有唯一解.') X=zeros(n,1); X(n)=b(n)/A(n,n);for k=n-1:-1:1 X(k)=(b(k)-sum(A(k,k+1:n)*X(k+1:n)/A(k,k); end else disp('请注意：因为RA=RB<n，所以此方程组有无穷多解.')endend例3.2.2 用解上三角形线性方程组的MATLAB程序解方程组.解
6、 在MATLAB工作窗口输入程序>>A=5 -1 2 3;0 -2 7 -4;0 0 6 5;0 0 0 3;b=20; -7; 4;6; RA,RB,n,X=shangsan(A,b)运行后输出结果请注意：因为RA=RB=n，所以此方程组有唯一解.RA = RB =4, 4,n = 4,X =2.4 -4.0 -1.0 2.03.3 高斯（Gauss）消元法和列主元消元法及其MATLAB程序3.3.1 高斯消元法及其MATLAB程序用高斯消元法解线性方程组的MATLAB程序function RA,RB,n,X=gaus(A,b)B=A b; n=length(b); RA=ran
7、k(A); RB=rank(B);zhica=RB-RA;if zhica>0,disp('请注意：因为RA=RB，所以此方程组无解.')returnendif RA=RB if RA=ndisp('请注意：因为RA=RB=n，所以此方程组有唯一解.') X=zeros(n,1); C=zeros(1,n+1); for p= 1:n-1for k=p+1:n m= B(k,p)/ B(p,p); B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);endend b=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)
8、/A(n,n); for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)/A(q,q); endelse disp('请注意：因为RA=RB<n，所以此方程组有无穷多解.')endend例3.3.2 用高斯消元法和MATLAB程序求解下面的非齐次线性方程组，并且用逆矩阵解方程组的方法验证.解 在MATLAB工作窗口输入程序>> A=1 -1 1 -3; 0 -1 -1 1;2 -2 -4 6;1 -2 -4 1; b=1;0; -1;-1; RA,RB,n,X =gaus (A,b)运行后输出结果请注意：因为RA=R
9、B=n，所以此方程组有唯一解.X = 0 -0.5000 0.5000 0RA = 4RB = 4n = 43.3.2 列主元消元法及其MATLAB程序用列主元消元法解线性方程组的MATLAB程序function RA,RB,n,X=liezhu(A,b)B=A b; n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA;if zhica>0,disp('请注意：因为RA=RB，所以此方程组无解.')returnendif RA=RB if RA=ndisp('请注意：因为RA=RB=n，所以此方程组有唯一解.')
10、 X=zeros(n,1); C=zeros(1,n+1); for p= 1:n-1Y,j=max(abs(B(p:n,p); C=B(p,:);B(p,:)= B(j+p-1,:); B(j+p-1,:)=C;for k=p+1:n m= B(k,p)/ B(p,p); B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);endend b=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n); for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)/A(q,q); endelse disp(&#
11、39;请注意：因为RA=RB<n，所以此方程组有无穷多解.')endend例3.3.3 用列主元消元法解线性方程组的MATLAB程序解方程组.解 在MATLAB工作窗口输入程序>> A=0 -1 -1 1;1 -1 1 -3;2 -2 -4 6;1 -2 -4 1; b=0;1;-1;-1; RA,RB,n,X=liezhu(A,b)运行后输出结果请注意：因为RA=RB=n，所以此方程组有唯一解.RA = 4，RB = 4，n = 4，X =0 -0.5 0.5 03.4 LU分解法及其MATLAB程序3.4.1判断矩阵LU分解的充要条件及其MATLAB程序判断矩阵能
12、否进行LU分解的MATLAB程序function hl=pdLUfj(A)n n =size(A); RA=rank(A); if RA=ndisp('请注意：因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的秩RA如下：'), RA,hl=det(A); returnendif RA=n for p=1:n,h(p)=det(A(1:p, 1:p);, endhl=h(1:n);for i=1:nif h(1,i)=0disp('请注意：因为A的r阶主子式等于零，所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：'),hl;RA,ret
13、urnendend if h(1,i)=0 disp('请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：')hl;RAendend例3.4.1 判断下列矩阵能否进行LU分解，并求矩阵的秩.（1）；（2）；（3）.解 （1）在MATLAB工作窗口输入程序>> A=1 2 3;1 12 7;4 5 6;hl=pdLUfj(A)运行后输出结果为请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：RA = 3， hl = 1 10 -48（2）在MATLAB工作窗口输入程
14、序 >> A=1 2 3;1 2 7;4 5 6;hl=pdLUfj(A)运行后输出结果为请注意：因为A的r阶主子式等于零，所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：RA = 3， hl =1 0 12（3）在MATLAB工作窗口输入程序 >> A=1 2 3;1 2 3;4 5 6;hl=pdLUfj(A)运行后输出结果为请注意：因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的秩RA如下RA = 2， hl = 03.4.2 直接LU分解法及其MATLAB程序将矩阵进行直接LU分解的MATLAB程序function hl=zhjLU(
15、A)n n =size(A); RA=rank(A); if RA=ndisp('请注意：因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的秩RA如下:'), RA,hl=det(A);returnendif RA=n for p=1:nh(p)=det(A(1:p, 1:p);endhl=h(1:n);for i=1:nif h(1,i)=0disp('请注意：因为A的r阶主子式等于零，所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：'), hl;RAreturnendend if h(1,i)=0 disp('请注意：因为A
16、的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：')for j=1:nU(1,j)=A(1,j);endfor k=2:nfor i=2:n for j=2:n L(1,1)=1;L(i,i)=1; if i>jL(1,1)=1;L(2,1)=A(2,1)/U(1,1); L(i,1)=A(i,1)/U(1,1);L(i,k)=(A(i,k)- L(i,1:k-1)*U(1:k-1,k)/U(k,k);elseU(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j);endendendendhl;RA,U,Lendend例3.
17、4.3 用矩阵进行直接LU分解的MATLAB程序分解矩阵.解 在MATLAB工作窗口输入程序>> A=1 0 2 0;0 1 0 1;1 2 4 3;0 1 0 3; hl=zhjLU(A)运行后输出结果L = 1 0 0 0 0 1 0 0 1 2 1 0 0 1 0 1 hl = 1 1 2 4 请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：RA = 4 U = 1 0 2 0 0 1 0 1 0 0 2 1 0 0 0 2 3.4.4 判断正定对称矩阵的方法及其MATLAB程序判断矩阵是否是正定对称矩阵的MATLAB程序
18、function hl=zddc(A)n n =size(A);for p=1:nh(p)=det(A(1:p, 1:p);endhl=h(1:n);zA=A'for i=1:n if h(1,i)<=0 disp('请注意：因为A的各阶顺序主子式hl不全大于零，所以A不是正定的.A的转置矩阵zA和各阶顺序主子式值hl依次如下：'), hl;zA,returnendendif h(1,i)>0disp('请注意：因为A的各阶顺序主子式hl都大于零，所以A是正定的.A的转置矩阵zA和各阶顺序主子式值hl依次如下：') hl;zAend例3.4.
19、5 判断下列矩阵是否是正定对称矩阵：（1）；(2) ； (3) ；（4）.解 （1）在MATLAB工作窗口输入程序>> A=0.1 2 3 4;-1 2 -3 4;11 21 13 41;5 7 8 9;hl=zddc(A)运行后输出结果请注意： A不是对称矩阵请注意：因为A的各阶顺序主子式hl不全大于零，所以A不是正定的.A的转置矩阵zA和各阶顺序主子式值hl依次如下：zA = 1/10 -1 11 5 2 2 21 7 3 -3 13 8 4 4 41 9 hl = 1/10 11/5 -1601/10 3696/5 因此,即不是正定矩阵，也不是对称矩阵.（2）在MATLAB工
20、作窗口输入程序>> A=1 -1 2 1;-1 3 0 -3;2 0 9 -6;1 -3 -6 19,hl=zddc(A)运行后输出结果A = 1 -1 2 1 -1 3 0 -3 2 0 9 -6 1 -3 -6 19 请注意： A是对称矩阵请注意：因为A的各阶顺序主子式hl都大于零，所以A是正定的.A的转置矩阵zA和各阶顺序主子式值hl依次如下：zA = 1 -1 2 1 -1 3 0 -3 2 0 9 -6 1 -3 -6 19 hl = 1 2 6 24 （3）在MATLAB工作窗口输入程序>> A=1/sqrt(2) -1/sqrt(2) 0 0; -1/sq
21、rt(2) 1/sqrt(2) 0 0; 0 0 1/sqrt(2) -1/sqrt(2); 0 0 -1/sqrt(2) 1/sqrt(2), hl=zddc(A)运行后输出结果A= 985/1393 -985/1393 0 0 -985/1393 985/1393 0 0 0 0 985/1393 -985/1393 0 0 -985/1393 985/1393 请注意： A是对称矩阵请注意：因为A的各阶顺序主子式hl不全大于零，所以A不是正定的.A的转置矩阵zA和各阶顺序主子式值hl依次如下：zA = 985/1393 -985/1393 0 0 -985/1393 985/1393 0
22、 0 0 0 985/1393 -985/1393 0 0 -985/1393 985/1393 hl = 985/1393 0 0 0 可见，不是正定矩阵，是半正定矩阵；因为= T 因此，是对称矩阵.（4）在MATLAB工作窗口输入程序 >> A=-2 1 1;1 -6 0;1 0 -4;hl=zddc(A)运行后输出结果 A = -2 1 1 1 -6 0 1 0 -4请注意： A是对称矩阵请注意：因为A的各阶顺序主子式hl不全大于零，所以A不是正定的.A的转置矩阵zA和各阶顺序主子式值hl依次如下：zA = -2 1 1 hl = -2 11 -38 1 -6 0 1 0 -
23、4可见不是正定矩阵，是负定矩阵；因为 = T 因此，是对称矩阵.3.5 求解线性方程组的LU方法及其MATLAB程序3.5.1 解线性方程组的楚列斯基（Cholesky）分解法及其MATLAB程序例3.5.1 先将矩阵进行楚列斯基分解,然后解矩阵方程，并用其他方法验证.解 在工作窗口输入>>A=1 -1 2 1;-1 3 0 -3; 2 0 9 -6;1 -3 -6 19;b1=1:2:7; b=b1' R=chol(A);C=A-R'*R,R1=inv(R);R2=R1' x=R1*R2*b,Rx=Ab-x运行后输出方程组的解和验证结果x = Rx = 1
24、.0e-014 * C = 1.0e-015 * -8.0000 -0.7105 0 0 0 0 0.3333 -0.0833 0 -0.4441 0 0 3.6667 0.2220 0 0 0 0 2.0000 0.1332 0 0 0 0 3.5.2 解线性方程组的直接LU分解法及其MATLAB程序例3.5.2 首先将矩阵直接进行LU分解，然后解矩阵方程，.解 (1) 首先将矩阵直接进行LU分解.在MATLAB工作窗口输入程序>> A=1 0 2 0;0 1 0 1;1 2 4 3;0 1 0 3;b=1;2;-1;5; hl=zhjLU(A),A-L*U运行后输出LU分解请注
25、意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：L = 1 0 0 0 0 1 0 0 1 2 1 00 1 0 1hl = 1 1 2 4RA = 4U = 1 0 2 0 0 1 0 1 0 0 2 1 0 0 0 2分解为一个单位下三角形矩阵和一个上三角形矩阵的积 .(2)在工作窗口输入>> U=1 0 2 0;0 1 0 1;0 0 2 1;0 0 0 2; L=1 0 0 0;0 1 0 0;1 2 1 0;0 1 0 1;b=1;2;-1;5;U1=inv(U); L1=inv(L); X=U1*L1*b,x=Ab运行后
26、输出方程组的解X = 8.50000000000000 0.50000000000000 -3.75000000000000 1.500000000000003.5.3 解线性方程组的选主元的LU方法及其MATLAB程序例3.5.3 先将矩阵进行LU分解，然后解矩阵方程 其中，.解 方法1 根据（3.28）式编写MATLAB程序，然后在工作窗口输入>> A=0.1 2 3 4;-1 2 -3 4;11 21 13 41;5 7 8 9; b=1;2;-1;5; L U P=LU(A), U1=inv(U); L1=inv(L); X=U1* L1*P*bP = 0 0 1 0 0
27、1 0 0 1 0 0 0 0 0 0 1X =-1.2013 3.3677 0.0536 -1.4440运行后输出结果L = 1.0000 0 0 0 -0.0909 1.0000 0 0 0.0091 0.4628 1.0000 0 0.4545 -0.6512 0.2436 1.0000U =11.0000 21.0000 13.0000 41.0000 0 3.9091 -1.8182 7.7273 0 0 3.7233 0.05120 0 0 -4.6171方法2 根据（3.29）式编写MATLAB程序，然后在工作窗口输入>> A=0.1 2 3 4;-1 2 -3 4;
28、11 21 13 41;5 7 8 9;b=1;2;-1;5; F U=LU(A), U1=inv(U); F1=inv(F); X=U1*F1*bU=11.0000 21.0000 13.0000 41.0000 0 3.9091 -1.8182 7.7273 0 0 3.7233 0.0512 0 0 0 -4.6171运行后输出结果F=0.0091 0.4628 1.0000 0 -0.0909 1.0000 0 0 1.0000 0 0 0 0.4545 -0.6512 0.2436 1.0000X =-1.2013 3.3677 0.0536 -1.4440用LU分解法解线性方程组的
29、MATLAB程序function RA,RB,n,X,Y=LUjfcz(A,b)n n =size(A);B=A b; RA=rank(A); RB=rank(B); for p=1:nh(p)=det(A(1:p, 1:p);endhl=h(1:n);for i=1:nif h(1,i)=0disp('请注意：因为A的r阶主子式等于零，所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：') hl;RAreturnendendif h(1,i)=0 disp('请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依
30、次如下：')X=zeros(n,1); Y=zeros(n,1); C=zeros(1,n);r=1:1;for p=1:n-1max1,j=max(abs(A(p:n,p); C=A(p,:); A(p,:)= A(j+P1,:); C= A(j+P1,:);g=r(p); r(p)= r(j+P1); r(j+P1)=g;for k=p+1:nH= A(k,p)/A(p,p); A(k,p) = H; A(k,p+1:n)=A(k,p+1:n)- H* A(p,p+1:n);endendY(1)=B(r(1);for k=2:nY(k)= B(r(k)- A(k,1:k-1)* Y
31、(1:k-1); endX(n)= Y(n)/ A(n,n);for i=n-1:-1:1 X(i)= (Y(i)- A(i, i+1:n) * X (i+1:n)/ A(i,i);endendRA,RB,n,X,Y;3.6 误差分析及其两种MATLAB程序3.6.1 用MATLAB软件作误差分析例3.6.2 解下列矩阵方程，并比较方程（1）和（2）有何区别，它们的解有何变化.其中解 (1) 矩阵方程的系数矩阵为7阶希尔伯特（Hilbert）矩阵，我们可以用下列命令计算阶希尔伯特矩阵 >>h=hilb(n) % 输出h为n阶Hilbert矩阵在MATLAB工作窗口输入程序 >
32、> A=hilb(7);b=1;2;2;2;2;2;2;X=Ab运行后输出的解为 X =（-35 504 -1260 -4200 20790 -27720 12012. （2）在MATLAB工作窗口输入程序>> B =0.001,zeros(1,6);zeros(6,1),zeros(6,6);A=(B+hilb(7); b=1;2;2;2;2;2;2;X=Ab 运行后输出方程的解为 X=（-33 465 -966 -5181 22409 -29015 12413. 在MATLAB工作窗口输入程序 >> X =-33, 465,-966,-5181,22409,-
33、29015,12413' X1 =-35,504,-1260,-4200,20790,-27720,12012' wu=X1'- X'运行后输出方程（1）和（2）的解的误差为 .方程（1）和（2）的系数矩阵的差为，常数向量相同，则的解的差为.的微小变化，引起的很大变化，即对的扰动是敏感的.3.6.2 求P 条件数和讨论解的性态的MATLAB程序求P条件数和讨论解的性态的MATLAB程序function Acp=zpjxpb(A)Acw = cond (A, inf);Ac1= cond (A,1);Ac2= cond (A,2);Acf = cond (A,&#
34、39;fro');dA=det(A);if (Acw>50)&(dA<0.1)disp('请注意：AX=b是病态的，A的条件数,1条件数,2条件数, 弗罗贝尼乌斯条件数和A的行列式的值依次如下：') Acp=Acw Ac1 Ac2 Acf dA'elsedisp(' AX=b是良态的，A的条件数,1条件数,2条件数,弗罗贝尼乌斯条件数和A的行列式的值依次如下：')Acp=Acw Ac1 Ac2 Acf dA'end例3.6.3 根据定理3.10，讨论线性方程组解的性态，并且求出的4种条件数.其中（1）为7阶希尔伯特矩阵
35、；（2）.解 （1）首先将求P条件数和讨论解的性态的MATLAB程序保存名为zpjxpb.m 的M文件，然后在MATLAB工作窗口输入程序>> Acp =zpjxpb(hilb(7); Acp',det(hilb(7)运行后输出结果请注意：AX=b是病态的，A的条件数,1条件数,2条件数, 弗罗贝尼乌斯条件数和A的行列式的值依次如下：ans = 1.0e+008 * 9.8519 9.8519 4.7537 4.8175 0.0000ans = 4.8358e-025（2）在MATLAB工作窗口输入程序 >> A=2 3 -1 5;3 1 2 -7;4 1 -3
36、 6;1 -2 4 -7;Acp=zpjxpb(A); Acp'运行后输出结果AX=b是良态的，A的条件数,1条件数,2条件数, 弗罗贝尼乌斯条件数和A的行列式的值依次如下：ans = 14.1713 19.4954 8.2085 11.4203 327.00003.6.3 用P范数讨论解和的性态的MATLAB程序用P范数讨论解和的性态的MATLAB程序function Acp=zpjwc(A,jA,b,jb,P)Acp = cond (A,P);dA=det(A); X=Ab;dertaA=A-jA; PndA=norm(dertaA, P);dertab=b-jb;Pndb=nor
37、m(dertab, P);if Pndb>0jX=Ajb; Pnb= norm(b, P);PnjX = norm(jX,P); dertaX=X-jX; PnjdX= norm(dertaX, P);jxX= PnjdX/PnjX; PnjX = norm(jX,P); PnX = norm(X,P); jxX= PnjdX/PnjX; xX= PnjdX/PnX; Pndb=norm(dertab,P); xAb=Pndb/Pnb;Pnbj=norm(jb,P); xAbj=Pndb/Pnbj; Xgxx= Acp*xAb;endif PndA>0jX=jAb; dertaX=
38、X-jX;PnX = norm(X,P); PnjdX= norm(dertaX, P); PnjX = norm(jX,P); jxX= PnjdX/PnjX;xX= PnjdX/PnX;PnjA=norm(jA,P); PnA=norm(A,P); PndA=norm(dertaA,P);xAbj= PndA/PnjA;xAb= PndA/PnA;Xgxx= Acp*xAb; endif (Acp >50)&(dA<0.1)disp('请注意：AX=b是病态的，A的P条件数Acp,A的行列式值dA，解X，近似解jX，解的相对误差jxX，解的相对误差估计值Xgxx
39、，b或A的相对误差xAb依次如下：') Acp,dA,X',jX',xX',jxX',Xgxx',xAb',xAbj'elsedisp('请注意： AX=b是良态的，A的P条件数Acp,A的行列式值dA，解X，近似解jX，解的相对误差jxX，解的相对误差估计值Xgxx，b或A的相对误差xAb依次如下：')Acp,dA,X',jX',xX',jxX',Xgxx',xAb',xAbj'end例3.6.4 根据定理3.10，讨论线性方程组解的性态，并利用（3.32）
40、式讨论当的每个元都取4位有效数字时，其解的相对误差.其中为7阶希尔伯特矩阵，.解 （1）取范数和条件数，线性方程组的不变时，取范数和条件数，系数矩阵为7阶希尔伯特矩阵，中的每个元素取4位有效数字.用P范数讨论解和的性态的MATLAB程序保存名为zpjwc.m的文件,然后在工作窗口输入MATLAB程序>> jA =1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769;A=hilb(7);b=1;1/3;4;2;2;2;2; jb=1;0.3333