数值计算方法课后实验

实验3.1Gauss消去法的数值稳定性实验

实验目的：观察和理解高斯消元过程中出现小主元（即la*%。很小）时引起方程组解数值不稳

定性.

实验内容：求解方程组Ax=b，其中

0.3X10-1559.143r59.17-

5.291-6.130--1246.78

(1)A,=9=;

11.29521

-1211..2.

■10-70I-■8

-32.099999999999625.900000000001

(2)A—瓦=5

25-15—1

0102_1

实验要求：

(1)计凫矩阵的条件数，判断系数矩阵是良态的还是病态的.

(2)用高斯列主元消去法求得L和U及解向量xl,x2.

(3)用不选主元的高斯消去法求得L和U及解向量xl,x2.

(4)观察小主元并分析对计算结果的影响.

解：⑴

cond(Al)=68.4296

cond(A2)=8.9939

Al矩阵条件数远大于1,A1矩阵为病态;A2矩阵条件数没有远远大于1,A2矩阵为良态。

(2)高斯列主元消去法：

程序如下

Clear

Al=(0.3*power(10,-15)Z59.14Z3Zl；5.291/-6.130r-l/2;11.2/9,5,2;1,2,1,1];

bl=(59.17;46.18;l;2];%Al=[lz2Z-1,1;1,1,2,-1;3,-1,1;2Zlf3Z-1];%

bl=[l;-2;6;-l];A2=[10,-7Z0,1；-3,2.099999999999,6,2;5,-1,5,-1;0,1,0,2];

b2=(8;5.900000000001;5;l];A=A1;b=bl;n=input(•0=');

fork=l:n-1[a,b3]=max(A(k:n,k));A([kk-l+b3],:)=A([k-l+b3k]z:);

b([kk-l+b3],:)=b([k-l+b3k]z:);ifA(k,k)〜=0

A(k+1:n,k)=A(k+l:n,k)./A(k,k);A(k+l:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:nzk

)*A(k,k+1:n);elsestopendendA;b;L=tril(Az-1);U=triu(A,0);

fori=l:1:nL(i,i)=l;endL;forj=l:n-1b(j)=b(j)/L(j,j);

b(j+l:n)=b(j+l:n)-b(j)*L(j+l:nzj);endb(n)=b(n)/L(nzn);b;y=b;

forj=n：-i：2y(i)=y(D/ud,i)；v(1：1-1)=v(1：1-1)-y(i)*u(i：i-1,j)；

endy(l)=y(l)/U(l,l);

结果如下：

方程组一：

LI:UI:

U=

L=

11.19009.00005.00002.0000

1.0000000

059.14203.00001.0000

0.00001.000000

000.49390.8012

0.08930.02021.00000

0005.8496

0.4734-0.1755-5.75281.0000

解向量Xl结果如下：

x=

3.8457

1.6020

-15.4761

10.3086

方程组二结果如下：

L2:U2:

u=

1.000000010.0000-7.000001.0000

0.50001.00000002.50005.0000-1.5000

-0.3000-0.00001.00000006.00002.3000

00.4000-0.33331.00000003.366；|

解向量x2结果如下：

x2=

0.0000

-1.0000

1.0000

1.0000

(3)不选主元的分解程序如下：

Clear

Al=[0.3*power(10,-15),59.14,3,1；5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1]；

bl=(59.17;46.18;l;2];%Al=[6,2,l,-l;2,4,l,0;l,l/4/-l；-l/0,-1,3];%bl=[6

;-l;5;-5];A2=[10,-7,0,1;-3,2.099999999999,6,2;5,-1,5,-1;0,1,0,2];b2=[

8;5.900000000001;5;l];A=A2;b=b2;n=input('n=');fork=l:n-l

A(k+1:nzk)=A(k+1:n,k)/A(kzk);A(k+1:nzk+1:n)=

A(k+l:n,k+l:n)-A(k+l:n,k)*A(k,k+l:n);endL=tril(A,-1);U=triu(Az0);

fori=l:1:nL(i,i)=1;endL;forj=l:n-1b(j)=b(j)/L(j,j);

b(j+l:n)=b(j+l:n)-b(j)*L(j+l:n,j);endb(n)=b(n)/L(nzn);b;y=b;

forj=n：-i：2y(i)=y(i)/u(i,i)；v(i：i-i)=v(i：i-i)-y('i)*u(i：i-i/i)；

endy(1)=y(1)/U(1,1);

方程组一结果:

LI:UI:

u=

0.00001.00000011.20009.00005.00002.0000

0.4724-0.17551.00000059.14003.00001.0000

1.000000000-2.83541.2307

0.08930.0202-0.17381.00000001.0151

解向量xl结果如卜.:

?3.679R

1.0004

0

0

方程组二结果：

L2:U2:

U=

1.0000000

-0.3000-0.00001.0000010.0000-7.000001.0000

0.50001.00000002.50005.0000-1.5000

00.4000-0.33331・0000006.00002.3000

0003.3667

解向量x2结果如下:

-0.0000

-1.0000

0.9999

1.0001

(4)观察方程在两种不同方法下的结果可知：

由于计算机字长是•定的，小主元会造成大数除以小数的结果超出字长，结果发生很大的变

化。

实验3.2方程组的性态和条件数实验

实验目的：

理解条件数的意义和方程组的性态对解向量的影响。

实验内容：

已知两个方程组八卜功和A2x=b,其中

ars=ans

N=8时，Al条件数为：.'A2条件数为：k526，e+1°

结论：随着n的增大，条件数也逐渐增大，矩阵病态更加严重；n=4,6,8时，矩

阵三种条件数非常大，都为病态矩阵。

(2)产5时，解向量XI与X2

1.0000

1.0000

1.0000

1.0000

1.0000

结果如下：A1结果为

1.0000

1.0000

1.0000

1.0000

1.0000

A2结果为:

(3)N=5时,对A1中a22项施加扰动leT4,解向量X为:

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

对A1中a66项施加扰动leT4,解向量X为:

xll=

1.0000

1.0000

1.0000

1.0000

1.0000

1.0000

对b元素施加扰动le-4,解向量X为：

0.7998

1.8458

-0.4249

2.1961

0.5000

1.0793

(4)取n=6不变,对A2中a22和a66施加扰动le-7,解向量分别为：

x2=x2=

1.00001.0000

1.00001.0000

1.00001.0000

1.00001.0000

1-00001.0000

1.000。i.oooo|

(5)对A1,A2的元素项进行微小扰动时，对解向量无明显影响；

对b元素进行微小扰动时，对解向量的影响较大；

结论：对Ax=b进行扰动时，解向量更容易受b元素的扰动影响，不容易受A元

素的扰动影响。

⑹求臼

1kL

N=5时，Al中a22和a66扰动对应的范数比分别为：

kl=k2=

2.7900e-101.3727e-10

实验3.3平方根法

实验目的：

掌握用平方根法、改进的平方根法求解系数矩阵正定的方程组的方法。

实验内容：

分别用平方根法、改进的平方根法求解方程组Ax=b,其中

「11/21/3…1/n-

1/21/31/4…l/(n+l)

A=1/31/41/5…l/(n+2)

■*••

■■••

[1/n1

/(n+1)l/(〃+2)•••1/(2〃-1).

nnT

b=，…，2a叼•

)—i-

实验要求：

(1)对于不同的n,编辑matlab程序求解方程组。

(2)与数学软件自带函数求解计完结果对比。

(3)与选主元的高斯消去法对比计算量。

解(1)对于不同的n,编辑matlab程序求解方程组。

程序如下：

平方根法程序

Clear

n=input('n=1);A=ones(n,n);

fori=l:1:n

forj=l:1:n

A(i,j)=(l/(i+j-l));

end

end

b=sum(A,2);A;b;%A=(4z-l/l;-l/4.25,2.75;lz2.75z3.5];%b=[6;-0.5;1.25];

fork=l:n

A(k,k)=sqrt(A(k,k));

A(kI1:nzk)-A(k11:n,k)/A(k,k);

forj=k+l:n

A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k)；

end

end

L=tril(A,0);

forj=l:n-l

b(j)=b(j)/L(j,j);

b(j+l:r))=b(j+l:n)-b(j)*L(j+l:n,j);

end

b(n)=b(n)/L(n,n);b;y=b;U=L*;

forj=n:-1:2

y(j)=y(j)/U(jzj)；y(l:j-l)=y(l:j-l)-y(j)*U(1:j-1,j);

end

y(l)=y(l)/U(l,l);

改进的平方根法程序：

Clear

A=ones(n,n);n=input(*n=*);

fori=l:1:n

for-j=l:1:n

A(i,j)=(l/(i+j-l));

end

end

b=Sum(A,2);%A=[4,-2,4,2;-2,10,-2,-7;4r-2,8f4;2,-7,4,7];

%b=(8;2;16;6];

forj=l:1:n

v=zeros(j-1,1);

fori=l:j-1

v(i)=A(j,i)*A(i,i);

end

A(j+l:nzj)=(A(j+l:nzj)-A(j+l:n,l:j-1)*v(l:j-1))/A(j,j);

end

A;L=tril(A,0);D=diag\A);

forj=l:n

L(jzj)=l；

end

L;forj=l:n-1

b(j)=b(j)/L(j,j);

b(j+l:n)=b(j+l:n)-b