MATLAB实验一 解线性方程组的直接法

实 验 报 告 课程名称 数值分析 实验项目 解线性方程组的直接法 专业班级 姓 名 学 号 指导教师 成 绩 日 期 月 日 一. 实验目的 1、掌握程序的录入和matlab的使用和操作； 2、了解影响线性方程组解的精度的因素——方法与问题的性态。 3、学会Matlab提供的“\\”的求解线性方程组。 二. 实验要求 1、按照题目要求完成实验内容； 2、写出相应的Matlab 程序； 3、给出实验结果(可以用表格展示实验结果)； 4、分析和讨论实验结果并提出可能的优化实验。 5、写出实验报告。 三. 实验步骤 1、用LU分解及列主元高斯消去法解线性方程组 ?7?10???32.099999a)?5?1??21?1??x1??8??????62??x2??5.900001???， 5?1??x3??5??????????02??x4??1?0输出Ax?b中系数A?LU分解的矩阵L和U，解向量x和det(A);用列主元法的行交换次序解向量x和求det(A)；比较两种方法所得结果。 2、用列主高斯消元法解线性方程组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、线性方程组Ax?b的A和b分别为 ?10??7A??8??7?77??32????565??23?, b????610933????31?5910????8则解x?(1,1,1,1,)T. 用MATLAB内部函数求det(A)和A的所有特征值和cond(A)2. 若令 78.17.2??10??5??7.085.046, A??A???85.989.899???6.99599.98???求解(A??A)(x??x)?b，输出向量?x和?x2，从理论结果和实际计算两方面分析线性方程组Ax?b解的相对误差?x4、 希尔伯特矩阵Hn?(hij)?Rn?n2/x2以及A的相对误差?A2/A2的关系。 1， i?j?1，其中hij?（1）分别对n?2,3,?,6计算cond(Hn)?，分析条件数作为n的函数如何变化。（2）令x?(1,1,?,1,)T?Rn，计算bn?Hnx，然后用高斯消去法解线性方程组Hnx?bn求出x，计算剩余向量rn?bn?Hnx以及?x?x?x。分析当n增加时解x分量的有效位数如何随n变化，它与条件数有何关系？当n多大时x连一位有效数字也

没有了？ 将每种情形的两个结果进行表格对比，如： n=6时： GAUSS列主消去法求得的x 四、实验结果 五、讨论分析 （对上述算例的计算结果进行比较分析，主要说清matlab的算符与消去法的适用范围不同，自己补充） 六、改进实验建议 （自己补充） x的有效数字 1.列主元的高斯消去法 利用列主元的高斯消去法matlab程序源代码： 首先建立一个gaussMethod.m的文件，用来实现列主元的消去方法。 function x=gaussMethod(A,b) %高斯列主元消去法，要求系数矩阵非奇异的， % n = size(A,1); if abs(det(A))<= 1e-8 error('系数矩阵是奇异的'); return; end % for k=1:n ak = max(abs(A(k:n,k))); index = find(A(:,k)==ak); if length(index) == 0 index = find(A(:,k)==-ak); end %交换列主元 temp = A(index,:);

A(index,:) = A(k,:); A(k,:) = temp; temp = b(index);b(index) = b(k); b(k) = temp; %消元过程 for i=k+1:n m=A(i,k)/A(k,k); %消除列元素 A(i,k+1:n)=A(i,k+1:n)-m*A(k,k+1:n); b(i)=b(i)-m*b(k); end end %回代过程 x(n)=b(n)/A(n,n); for k=n-1:-1:1; x(k)=(b(k)-A(k,k+1:n)*x(k+1:n)')/A(k,k); end x=x'; end 然后调用gaussMethod函数，来实现列主元的高斯消去法。在命令框中输入下列命令： 输出结果如下： 利用LU分解法及matlab程序源代码： function [L,U]=myLU(A) %实现对矩阵A的LU分解，L为下三角矩阵 A[n,n]=size(A);

L=zeros(n,n); U=zeros(n,n); for i=1:n L(i,i)=1; end for k=1:n for j=k:n U(k,j)=A(k,j)-sum(L(k,1:k-1).*U(1:k-1,j)'); end for i=k+1:n L(i,k)=(A(i,k)-sum(L(i,1:k-1).*U(1:k-1,k)'))/U(k,k); end end 在命令框中输入下列命令： >> a=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2] a = 10.0000 -7.0000 0 1.0000 -3.0000 2.1000 6.0000 2.0000 5.0000 -1.0000 5.0000 -1.0000 2.0000 1.0000 0 2.0000 >> [l,u]=lu(a) l = 1.0000 0 0 0 -0.3000 -0.0000 1.0000 0 0.5000 1.0000 0 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 >> b=[8 5.900001 5 1]' b = 8.0000 5.9000 5.0000 1.0000 >> y=l\\b y = 8.0000 1.0000 8.3000

5.0800 >> x1=U\\x x1 = 0.0000 -1.0000 1.0000 1.0000 >>det1= det(a) det1 =-762.0001 2、（1）在MATLAB窗口： >>A=[3.01 6.03 1.99;1.27 4.16 -1.23;0.987 -4.81 9.34] A = 3.0100 6.0300 1.9900 1.2700 4.1600 -1.2300 0.9870 -4.8100 9.3400 >> b=[1 1 1]' b = 1 1 1 >>[x1,det1,index]=Gauss (A,b) x1 = 1.0e+03 * 1592.599624841381 -631.9113762025488 -493.6177247593899 det1 = -0.0305 index = 1

(2) 在MATLAB窗口： >> A=[3.00 6.03 1.99;1.27 4.16 -1.23;0.990 -4.81 9.34] A = 3.0000 6.0300 1.9900 1.2700 4.1600 -1.2300 0.9900 -4.8100 9.3400 >> b=[1 1 1]' b = 1 1 1 >>[x2,det2,index]=Gauss5555(A,b) x2 = 119.5273 -47.1426 -36.8403 det2 = -0.4070 index = 1 3、在MATLAB窗口： A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10]; b=[32 23 33 31]'; x=A\\b b1=[32.1 22.9 33.1 30.9]'; x1=A\\b1 A1=[10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98]; x2=A1\\b delta_b=norm(b-b1)/norm(b) delta_A=norm(A-A1)/norm(A) delta_x1=norm(x-x1)/norm(x) delta_x2=norm(x-x2)/norm(x)

cond_A=cond(A) x = 1.0000 1.0000 1.0000 1.0000 x1 = 9.2000 -12.6000 4.5000 -1.1000 x2 = -9.5863 18.3741 -3.2258 3.5240 delta_b = 0.0033 delta_A = 0.0076 delta_x1 = 8.1985 delta_x2 = 10.4661 cond_A = 2.9841e+03 3、在MATLAB窗口： A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10]; b=[32 23 33 31]'; x=A\\b b1=[32.1 22.9 33.1 30.9]'; x1=A\\b1 A1=[10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98]; x2=A1\\b delta_b=norm(b-b1)/norm(b)

delta_A=norm(A-A1)/norm(A) delta_x1=norm(x-x1)/norm(x) delta_x2=norm(x-x2)/norm(x) cond_A=cond(A) x = 1.0000 1.0000 1.0000 1.0000 x1 = 9.2000 -12.6000 4.5000 -1.1000 x2 = -9.5863 18.3741 -3.2258 3.5240 delta_b = 0.0033 delta_A = 0.0076 delta_x1 = 8.1985 delta_x2 = 10.4661 cond_A = 2.9841e+03 4、 k=13 for n=2:6 a=hilb(n); co(n)=cond(a,inf); end x=1:6; plot(x,co);

b=zeros(k); x11=b; x0=b; r=b; for i=2:k x=(linspace(1,1,i))'; x0(1:i,(i-1))=x; H=hilb(i); b0=H*x; b(1:i,(i-1))=b0; x1=gauss2(H,b0); r(1:i,(i-1))=b0-H*x1; x11(1:i,(i-1))=x1; end dx=x11-x0; 结果如下： co=[0,27, 748, 28375,943656, 29070279]可见，条件数随着n的增大，急剧增加，如图1所示。 将（2）求得的结果（dx，x11）整理得 n=2时： x 1 1 n=3时： Δx 4.44E-16 -6.66E-16 rn 0 0 有效数字 16 15 x 1 1 1 n=4时： Δx -1.33E-15 9.55E-15 -9.99E-15 rn -2.22E-16 0 0 有效数字 15 14 14 x 1 1 1 1 n=5时： Δx -2.35E-14 2.56E-13 -6.11E-13 3.96E-13 rn 0 0 1.11E-16 0 有效数字 14 13 12 13 x Δx rn 有效数字

1 1 1 1 1 n=6时： -1.21E-14 6.97E-14 1.18E-13 -6.17E-13 4.59E-13 4.44E-16 0 2.22E-16 0 -1.11E-16 14 13 13 12 13 x 1 1 1 1 1 1 n=7时： Δx -9.28E-13 2.67E-11 -1.82E-10 4.75E-10 -5.26E-10 2.08E-10 rn 0 0 0 0 0 0 有效数字 12 11 10 10 9 10 x 1 1 1 1.00000001 0.99999997 1.00000002 0.99999999 n=8时： Δx -9.26E-12 3.71E-10 -3.59E-09 1.40E-08 -2.57E-08 2.22E-08 -7.30E-09 rn 0 0 2.22E-16 0 0 1.11E-16 0 有效数字 11 10 9 8 8 8 8 x 1 1 0.99999998 1.00000011 0.9999997 1.00000042 0.9999997 1.00000008 n=9时： Δx -2.82E-11 1.53E-09 -2.01E-08 1.09E-07 -2.95E-07 4.19E-07 -2.99E-07 8.44E-08 rn 4.44E-16 0 2.22E-16 0 -2.22E-16 -1.11E-16 1.11E-16 0 有效数字 11 9 8 7 7 7 7 7

x 1 1.00000002 0.99999968 1.00000228 0.99999164 1.00001706 0.99998042 1.00001182 0.99999708 n=10时： Δx -2.75E-10 1.90E-08 -3.20E-07 2.28E-06 -8.36E-06 1.71E-05 -1.96E-05 1.18E-05 -2.92E-06 rn 4.44E-16 0 0 4.44E-16 0 0 -1.11E-16 0 0 有效数字 10 8 7 6 5 5 5 5 6 x 1 1.00000008 0.99999835 1.00001491 0.99992911 1.00019443 0.99968164 1.00030715 0.99983898 1.00003537 n=11时： Δx -9.05E-10 7.76E-08 -1.65E-06 1.49E-05 -7.09E-05 0.000194426 -0.000318365 0.000307149 -0.000161019 3.54E-05 rn 4.44E-16 0 0 4.44E-16 0 0 -1.11E-16 0 0 0 有效数字 9 7 6 5 4 4 4 4 4 5 x 0.999999995 1.000000539 0.999986041 1.000155875 0.999072325 1.003258718 0.992910161 1.009658807 0.991981908 1.003707654 0.999267974 Δx -5.16E-09 5.39E-07 -1.40E-05 0.00015588 -0.0009277 0.00325872 -0.0070898 0.00965881 -0.0080181 0.00370765 -0.000732 rn 0 -4.44E-16 0 0 0 0 0 -2.22E-16 0 0 -2.22E-16 有效数字 8 6 5 4 3 3 2 2 3 3 3

n=12时： x 0.999999965 1.000004409 0.999861744 1.001879711 0.986241983 1.060375974 0.831939173 1.303973363 0.643862006 1.260684018 0.891666924 1.019510753 n=13时： Δx -3.49E-08 4.41E-06 -0.0001383 0.00187971 -0.013758 0.06037597 -0.1680608 0.30397336 -0.356138 0.26068402 -0.1083331 0.01951075 rn 8.88E-16 4.44E-16 0 0 0 2.22E-16 2.22E-16 0 1.11E-16 2.22E-16 1.11E-16 0 有效数字 8 6 4 3 2 2 1 1 1 1 1 2 x 1.000000069 0.999989224 1.000413538 0.993147495 1.061246208 0.669261244 2.149074131 -1.65417119 5.11854808 -3.243049939 3.783004896 -0.051792817 1.174329117 Δx 6.89E-08 -1.08E-05 0.00041354 -0.0068525 0.06124621 -0.3307388 1.14907413 -2.6541712 4.11854808 -4.2430499 2.7830049 -1.0517928 0.17432912 rn 8.88E-16 -4.44E-16 -2.22E-16 2.22E-16 -2.22E-16 2.22E-16 0 0 -1.11E-16 1.11E-16 0 1.11E-16 1.11E-16 有效数字 7 5 4 3 2 1 0 0 0 0 0 0 1 五、分析讨论： 实验的数学原理很容易理解，也容易上手。把运算的结果带入原方程组，可以发现符合的还是比较好。这说明列主元消去法计算这类方程的有效性。当A可逆时，能够将计算进行到底，列主元法就能确保算法的稳定，而且计算量不大。直接三角消去过程，实质上

是将A分解为两个三角矩阵的乘积A=LU，并求解Ly=b的过程。回带过程就是求解上三角方程组Ux=y。所以在实际的运算中，矩阵L和U可以直接计算出，而不需要任何中间步骤，从而在计算过程中将高斯消去法的步骤进行了进一步的简略，大大提高了运算速度，这就是三角分解法。 通过以上的计算比较，2.题方程组具有严重的病态性。当系数矩阵有微小的变化时，wucha =-401.8918 159.5435 124.6330，所得的解与原方程组的解有很大的相对误差。1题方程组中当系数矩阵A和b有微小变化时，wucha =0 0 0 0，所得的解与方程组的解没有相对误差。所以1题方程组是良性的。用MATLAB内部函数inv通过求逆矩阵，然后通过x=inv（A）*b也可以求出方程组的解，但是没有列主元高斯消去法具有良好的稳定性。det函数求方程组系数矩阵的行列式时所得结果和高斯消去法和三角法所得结果相同，具有方便快捷的优点。题四可以看出，条件数越大，有效位数越少，当n=13时，出现有效位数为0的情况。

附：高斯列主消去法源程序代码

function [x,det,index]=Gauss(A,b) % ?ó??D?·?3ì×éμ?áD?÷?aGauss??è￥·¨￡????D % A --- 方程组矩阵 % b --- 方程组右端 % x --- 方程组的解 % det ---方程组行列式

% index --- index=0表示求解失败，index=1表示求解成功|

[n,m]=size(A); nb=length(b); if n~=m

error('The rows and columns of matrix A must be equal!'); return; end if m~=nb

error('The columns of A must be equal the dimension of b!'); return; end

index=1; det=1; x=zeros(n,1); for k=1:n-1 % 选主元 a_max=0; for i=k:n

if abs(A(i,k))>a_max a_max=abs(A(i,k)); r=i; end end

if a_max<1e-20 index=0; return; end % 交换两行 if r>k for j=k:n

z=A(k,j); A(k,j)=A(r,j); A(r,j)=z; end

z=b(k); b(k)=b(r); b(r)=z; det=-det; end

% 消元过程 for i=k+1:n

m=A(i,k)/A(k,k); for j=k+1:n

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

b(i)=b(i)-m*b(k); end

det=det*A(k,k); end

det=det*A(n,n); % 回代过程

if abs(A(n,n))<1e-20 index=0; return; end

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

b(k)=b(k)-A(k,j)*x(j); end

x(k)=b(k)/A(k,k); end

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

b(i)=b(i)-m*b(k); end

det=det*A(k,k); end

det=det*A(n,n); % 回代过程

if abs(A(n,n))<1e-20 index=0; return; end

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

b(k)=b(k)-A(k,j)*x(j); end

x(k)=b(k)/A(k,k); end

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