数值分析第07次作业

实验七 直接解线性方程组的数值实验

实验1：

题目：用列主元Gauss消元法求解以下方程组：

?0.7290.810.9??x1??0.6867???????111x?0.83382???????1.3311.211.1??x??1.000????3???

编程及结果：

编程：

fprintf('增广矩阵A')

A=[0.729 0.81 0.9 0.6867;1 1 1 0.8338;1.331 1.21 1.1 1.000] tempo=A(3,:);A(3,:)=A(1,:);A(1,:)=tempo; fprintf('第一次选主元后的矩阵') A

A(2,:)=A(2,:)-A(1,:)*A(2,1)/A(1,1); A(3,:)=A(3,:)-A(1,:)*A(3,1)/A(1,1); fprintf('第一次消元后的矩阵') A

tempo=A(3,:);A(3,:)=A(2,:);A(2,:)=tempo; fprintf('第二次选主元后的矩阵') A

A(3,:)=A(3,:)-A(2,:)*A(3,2)/A(2,2); fprintf('第二次消元后的矩阵') A

%回代求解方程组 x(3)=A(3,4)/A(3,3);

x(2)=(A(2,4)-A(2,3)*x(3))/A(2,2);

x(1)=(A(1,4)-A(1,3)*x(3)-A(1,2)*x(2))/A(1,1); x

结果： 增广矩阵A A =

Columns 1 through 2

0.7290 0.8100 1.0000 1.0000 1.3310 1.2100

Columns 3 through 4

0.9000 0.6867 1.0000 0.8338 1.1000 1.0000

第一次选主元后的矩阵 A =

Columns 1 through 2

1.3310 1.2100 1.0000 1.0000 0.7290 0.8100

Columns 3 through 4

1.1000 1.0000 1.0000 0.8338 0.9000 0.6867

第一次消元后的矩阵 A =

Columns 1 through 2

1.3310 1.2100 0 0.0909 0 0.1473

Columns 3 through 4

1.1000 1.0000 0.1736 0.0825 0.2975 0.1390

第二次选主元后的矩阵 A =

Columns 1 through 2

1.3310 1.2100 0 0.1473 0 0.0909

Columns 3 through 4

1.1000 1.0000 0.2975 0.1390 0.1736 0.0825

第二次消元后的矩阵 A =

Columns 1 through 2

1.3310 1.2100 0 0.1473 0 0

Columns 3 through 4

1.1000 1.0000 0.2975 0.1390 -0.0101 -0.0033 x =

Columns 1 through 2

0.2245 0.2814

Column 3

0.3279

实验2：

题目：写用追赶法解三对角方程组的程序，并解下列方程组：

?2x1?x2?5??x?2x?x??12?123 ???x2?2x3?x4?11???x3?2x4??1

编程及结果：

编程：

A=[2,-1,0,0;-1,2,-1,0;0,-1,2,-1;0,0,-1,2]; b=[5,-12,11,-1]';

[l,u,p]=lu(A);%对矩阵A进行LU分解,l单位下三角，u上三角，p是置换矩阵 %l,u

y(1)=b(1);

y(2)=b(2)-l(2,1)*y(1);

y(3)=b(3)-l(3,1)*y(1)-l(3,2)*y(2);

y(4)=b(4)-l(4,1)*y(1)-l(4,2)*y(2)-l(4,3)*y(3); x(4)=y(4)/u(4,4);

x(3)=(y(3)-u(3,4)*x(4))/u(3,3);

x(2)=(y(2)-u(2,4)*x(4)-u(2,3)*x(3))/u(2,2);

x(1)=(y(1)-u(1,4)*x(4)-u(1,3)*x(3)-u(1,2)*x(2))/u(1,1); x

结果： x =

1.0000 -3.0000 5.0000 2.0000

实验3：

题目：用LU分解及列主元高斯消去法解线性方程组

?701??x1??8?10????x????32.099999625.900001??2???? ??5??15?1??x3??5??????x21021???4???输出Ax?b中系数A?LU分解的矩阵L和U，解向量x及det(A)；列主元法的行交换次序，解向量x及det(A)；比较两种方法所得的结果。

编程及结果：

1)LU分解

clc clear

A=[10 -7 0 1

-3 2.099999 6 2 5 -1 5 -1 2 1 0 2 ];

b=[8 5.900001 5 1];

[l,u,p]=lu(A);%对矩阵A进行LU分解,l单位下三角，u上三角，p是置换矩阵 l,u

y(1)=b(1);

y(2)=b(2)-l(2,1)*y(1);

y(3)=b(3)-l(3,1)*y(1)-l(3,2)*y(2);

y(4)=b(4)-l(4,1)*y(1)-l(4,2)*y(2)-l(4,3)*y(3); x(4)=y(4)/u(4,4);

x(3)=(y(3)-u(3,4)*x(4))/u(3,3);

x(2)=(y(2)-u(2,4)*x(4)-u(2,3)*x(3))/u(2,2);

x(1)=(y(1)-u(1,4)*x(4)-u(1,3)*x(3)-u(1,2)*x(2))/u(1,1); x

t=det(A) 结果： l =

1.0000 0 0 0 0.5000 1.0000 0 0 -0.3000 -0.0000 1.0000 0 0.2000 0.9600 -0.8000 1.0000 u =

10.0000 -7.0000 0 1.0000 0 2.5000 5.0000 -1.5000 0 0 6.0000 2.3000 0 0 0 5.0800

x = 0.1949 -0.7661 0.9695 0.6882 t =

-762.0001

2）列主元高斯消去法 clear clc

A=[10 -7 0 1 ; -3 2.099999 6 2 ; 5 -1 5 -1 ; b=[8 5.900001 5 1]; n=length(b); X=zeros(n,1); c=zeros(1,n); d1=0; t=det(A)

2 1 0 2 ];

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 A 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 >> t =

762.0001 A =

10.0000 -7.0000 0 1.0000 0 -0.0000 6.0000 2.3000 0 2.5000 5.0000 -1.5000 0 2.4000 0 1.8000 A =

10.0000 -7.0000 0 1.0000 0 2.5000 5.0000 -1.5000 0 0 6.0000 2.3000 0 0 -4.8000 3.2400 A =

10.0000 -7.0000 0 1.0000 0 2.5000 5.0000 -1.5000 0 0 6.0000 2.3000 0 0 0 5.0800 X =

0.0000 -1.0000 1.0000 1.0000 t =

762.0001

实验4：

题目：用列主元高斯消去法解线性方程组Ax?b：

?3.016.031.99??x1??1???????4.16?1.23??x2???1?; （1）?1.27?0.987?4.819.34??x??1????3????3.006.031.99??x1??1???????4.16?1.23??x2???1?; （2）?1.27?0.990?4.819.34??x??1????3???分别输出A,b,det(A)，解向量x，(1)中A的条件数，分析比较(1),(2)的计算结果。

编程及结果：

?3.016.031.99??x1??1???????4.16?1.23??x2???1?; （1）?1.27?0.987?4.819.34??x??1????3???编写程序： clear clc

A=[3.01 6.03 1.99; 1.27 4.16 -1.23; 0.987 -4.81 9.34] b=[1 1 1] n=length(b); X=zeros(n,1); c=zeros(1,n); d1=0;

fprintf('输出det(A)') det(A)

fprintf('输出条件数：1-条件数，2-条件数，无穷-条件数') [cond(A,1),cond(A),cond(A,inf)] 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

运行结果： A =

3.0100 6.0300 1.9900 1.2700 4.1600 -1.2300 0.9870 -4.8100 9.3400 b =

1 1 1

输出det(A) ans =

-0.0305

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

1.0e+004 *

5.5228 3.0697 5.6751

X =

1.0e+003 *

1.5926 -0.6319 -0.4936

?3.006.031.99??x1??1???????4.16?1.23??x2???1?; （2）?1.27?0.990?4.819.34??x??1????3???

编写程序： clear clc

A=[3.00 6.03 1.99; 1.27 4.16 -1.23; 0.990 -481 9.34] b=[1 1 1] n=length(b); X=zeros(n,1); c=zeros(1,n); d1=0;

fprintf('输出det(A)') det(A)

fprintf('输出条件数：1-条件数，2-条件数，无穷-条件数') [cond(A,1),cond(A),cond(A,inf)] 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

运行结果： 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

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

本文来源：https://www.bwwdw.com/article/ysfp.html