数值分析第07次作业(3)

end

for k=i+1:n for j=i+1:n

A(k,j)=A(k,j)-A(i,j)*A(k,i)/A(i,i); end

b(k)=b(k)-b(i)*A(k,i)/A(i,i); A(k,i)=0; end end

X(n)=b(n)/A(n,n); for i=n-1:-1:1 sum=0; for j=i+1:n

sum=sum+A(i,j)*X(j); end

X(i)=(b(i)-sum)/A(i,i); end X

运行结果： A =

3.0000 6.0300 1.9900 1.2700 4.1600 -1.2300 0.9900 -481.0000 9.3400 b =

1 1 1

输出det(A) ans =

-2.9610e+003

输出条件数：1-条件数，2-条件数，无穷-条件数 ans =

412.8103 289.3865 343.2779 X =

0.5343 -0.0065 -0.2833

分别输出A,b,det(A)，解向量x，(1)中A的条件数，分析比较(1),(2)的计算结果。

提高题:

题目：用高斯消元法、列主元的高斯消元法和追赶法分别求解三对角线性方程组：

?4?1??x1??7????x????14?15???2??????x3???13??14?1??????x?14?12???4??????x5??6??14?1????????14?1???x6???12????x??14??14?1???7????14?1???x8???4????x??5??14?1???9?????????14??x10???5???

Tx*?(2,1,?3,0,1,?2,3,0,1,?1)其中精确解为

编程及结果：

1) 高斯消元法 编写程序： clear clc

A=[ 4 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0

0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 4 ]; b=[7 5 -13 2 6 -12 14 -4 5 -5]; n=length(b); x=b;

for i=1:n-1 for j=i+1:n

mi=A(j,i)/A(i,i); b(j)=b(j)-mi*b(i); for k=i:n

A(j,k)=A(j,k)-mi*A(i,k); end [A b']; end end

x(n)=b(n)/A(n,n); for i=n-1:-1:1 s=0;

for j=i+1:n

s=s+A(i,j)*x(j); end

x(i)=(b(i)-s)/A(i,i); end x

运行结果： x =

Columns 1 through 6

2.0000 1.0000 -3.0000 0.0000

Columns 7 through 10

3.0000 -0.0000 1.0000 -1.0000

2) 列主元的高斯消元法 编写程序： clear clc

-2.0000 1.0000 A=[ 4 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 4 ]; b=[7 5 -13 2 6 -12 14 -4 5 -5]; n=length(b); X=zeros(n,1); c=zeros(1,n); d1=0;

for i=1:n-1

max=abs(A(i,i)); m=i;

for j=i+1:n

if max

for k=i:n

c(k)=A(i,k); A(i,k)=A(m,k); A(m,k)=c(k); end d1=b(i); b(i)=b(m); b(m)=d1; end

for k=i+1:n for j=i+1:n

A(k,j)=A(k,j)-A(i,j)*A(k,i)/A(i,i); end

b(k)=b(k)-b(i)*A(k,i)/A(i,i); A(k,i)=0; end end

X(n)=b(n)/A(n,n);

for i=n-1:-1:1 sum=0; for j=i+1:n

sum=sum+A(i,j)*X(j); end

X(i)=(b(i)-sum)/A(i,i); end X

运行结果： >> X =

2.0000 1.0000 -3.0000 0.0000 1.0000 -2.0000 3.0000 -0.0000 1.0000 -1.0000

3) 追赶法 编写程序： clear clc

A=[ 4 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 4 ];

d=[7 5 -13 2 6 -12 14 -4 5 -5];%定义三对角矩阵A的各组成单元。方程为Ax=d [l,u,p]=lu(A);

a=[-1 -1 -1 -1 -1 -1 -1 -1 -1 -1];% a为-1对角线元素(2~n)， b=[4 4 4 4 4 4 4 4 4 4 ];% b为A的对角线元素(1~n)

c=[-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ];% c为+1对角线元素(1~n-1)。

n=length(b);

u0=0;y0=0;a(1)=-1; L(1)=b(1)-a(1)*u0;

y(1)=(d(1)-y0*a(1))/L(1); u(1)=c(1)/L(1); for i=2:(n-1)

L(i)=b(i)-a(i)*u(i-1);

y(i)=(d(i)-y(i-1)*a(i))/L(i); u(i)=c(i)/L(i); end

L(n)=b(n)-a(n)*u(n-1);

y(n)=(d(n)-y(n-1)*a(n))/L(n); x(n)=y(n); for i=(n-1):-1:1

x(i)=y(i)-u(i)*x(i+1); end x

运行结果： >> x =

Columns 1 through 6

2.0000 1.0000 -3.0000

Columns 7 through 10

3.0000 -0.0000 1.0000

0.0000 1.0000 -1.0000 -2.0000