大连理工大学矩阵与数值分析上机作业代码

《矩阵与数值分析》上机作业实验报告
陈晓 学号：21204047 机械工程学院 教学班号：02 任课教师：金光日
习题 1
给定 n 阶方程组 Ax = b ，其中
⎛6 1 ⎞ ⎛7⎞ ⎜ ⎟ ⎜ ⎟ ⎜8 6 1 ⎟ ⎜ 15 ⎟ A = ⎜ ⋱ ⋱ ⋱ ⎟ ，b = ⎜ ⋮ ⎟ ⎜ ⎟ ⎜ ⎟ 8 6 1⎟ ⎜ ⎜ 15 ⎟ ⎜ ⎟ ⎜14 ⎟ 8 6⎠ ⎝ ⎝ ⎠
则方程组有解 x = (1,1,⋯ ,1) 。 对 n = 10 和 n = 84 ， 分别用 Gauss 消去法和列主元消去法解
T
方程组，并比较计算结果。
1.1 程序
（1）高斯消元法 n=10 时， >> [A,b]=CreateMatrix(10) A= 6 8 0 0 0 0 0 0 0 0 1 6 8 0 0 0 0 0 0 0 0 1 6 8 0 0 0 0 0 0 0 0 1 6 8 0 0 0 0 0 0 0 0 1 6 8 0 0 0 0 0 0 0 0 1 6 8 0 0 0 0 0 0 0 0 1 6 8 0 0 0 0 0 0 0 0 1 6 8 0 0 0 0 0 0 0 0 1 6 8 0 0 0 0 0 0 0 0 1 6
U=p{i}*U; for j=i+1:n l{i}(j,i)=-U(j,i)/Max; end U=l{i}*U; end %计算L L=eye(n); for i=1:n-2 for j=i+1:n-1 l{i}=p{j}*l{i}/p{j}; end L=l{i}*L; end L=inv(L*l{n-1}); %计算P P=eye(n); for i=1:n-1 P=p{i}*P; end %计算x x=U\(L\P*b);
-0.0209 0.0419 -0.0836 0.1665 -0.3303 0.6501 -1.2582 2.3487 -4.0263 5.3684 （2）列主元消去法 n=10 时， >> [A,b]=CreateMatrix(10) >> x=MainElementMethod(A,b) x= 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 n=84 时， >> [A,b]=CreateMatrix(84) >> x=MainElementMethod(A,b) x= 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
0.3111 -0.0882 0.3104 -0.0957 0.3313 -0.0472 n= 9 （5）取 x
(0)
=（1,1， ……，1,1）T， b =（1,1， ……，1,1）T,迭代误差取默认值
10
−5
，迭代次
数上限取默认值 50，使用 Jacobi 法进行迭代。 >> test2 >> b=zeros(20,1) >> x0=zeros(20,1) >> [x,n]=JacobiMethod(A,b,x0) x= 0.2766 0.2327 0.2159 0.2223 0.2228 0.2221 0.2222 0.2222 0.2222 0.2222 0.2222 0.2222 0.2222 0.2222 0.2221 0.2228 0.2223 0.2159 0.2327 0.2766 n= 19 （6）取 x
将 A 的主对角线元素成倍增长若干次，非主对角元素不变，每次用 Jacobi 法计算，要求迭 代误差满足 x ( k +1) − x ( k )
∞
< 10−6 ，比较收敛速度，分析现象并得出你的结论。
2.1 程序
(0) −5
（1）取 x
Байду номын сангаас
= 0 ， b =（1,1， ……，1,1）T,迭代误差取默认值
n= 9 （3）取 x
(0)
= 0 ， b =（1,0,1,0,……，1,0）T,迭代误差取默认值
10
−5
，迭代次数上限取默认
值 50，使用 Jacobi 法进行迭代。 >> test2 >> b=[1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0]' >> x0=zeros(20,1) >> [x,n]=JacobiMethod(A,b,x0) x= 0.3238 -0.0985 0.3116 -0.0881 0.3110
10
，迭代次数上限取默认值
50，使用 Jacobi 法进行迭代。 >> test2 >> b=ones(20,1) >> x0=zeros(20,1) >> [x,n]=JacobiMethod(A,b,x0) x= 0.2766 0.2327 0.2159 0.2223 0.2227 0.2221 0.2222 0.2222 0.2222 0.2222 0.2222 0.2222 0.2222 0.2222 0.2221 0.2227 0.2223 0.2159 0.2327 0.2766
习题 2
设有方程组 Ax = b ，其中 A ∈ R
20× 20
，
−0.5 −0.25 ⎛ 3 ⎜ 3 −0.5 ⎜ −0.5 ⎜ −0.25 −0.5 ⋱ ⎜ A=⎜ ⋱ ⋱ ⎜ ⋱ ⎜ ⎜ ⎜ ⎝
(a) 选 取 不 同 的 初 始 向 量 x
(0)
⋱ ⋱ ⋱ ⋱ ⋱
⎞ ⎟ ⎟ ⎟ ⋱ ⎟ ⋱ ⎟ ⋱ −0.5 −0.25 ⎟ ⎟ −0.5 3 −0.5 ⎟ −0.25 −0.5 3 ⎟ ⎠
10
−5
，迭代次数上限取默认
值 50，使用 Gauss-Seidel 法进行迭代。 >> test2 >> b=[1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0]' >> x0=zeros(20,1) >> [x,n]=Gauss_SeidelMethod(A,b,x0) x= 0.3238 -0.0985 0.3115 -0.0881 0.3110 -0.0889 0.3111 -0.0889 0.3111 -0.0889 0.3111 -0.0889 0.3111 -0.0889
1.3.2 GaussMethod.m function x=GaussMethod(A,b) %高斯消元法返回方程Ax=b的解
%x是方程的解 %A为方程的系数矩阵 %b为方程的常数矩阵 n=length(b); %定义单位下三角矩阵元胞数组，并将全部数组元素赋值为n阶单位阵 l=cell(1,n-1); l(:)={eye(n)}; U=A; for i=1:n-1 for j=i+1:n l{i}(j,i)=-U(j,i)/U(i,i); end U=l{i}*U; end %计算L和U L=inv(l{1}); for i=2:n-1 L=L/l{i}; end x=U\(L\b); .m 1.3.3 MainElementMethod MainElementMethod.m function x=MainElementMethod(A,b) %列主元消去法返回方程Ax=b的解 %x是方程的解 %A为方程的系数矩阵 %b为方程的常数矩阵 n=length(b); %定义单位下三角矩阵元胞数组，并将全部数组元素赋值为n阶单位阵 l=cell(1,n-1); l(:)={eye(n)}; p=cell(1,n-1); p(:)={eye(n)}; U=A; for i=1:n-1 %找出列主元所在行，并交换首行，得到Pi [Max,row]=max(U(i:n,i)); row=row+i-1; mid=p{i}(row,:); p{i}(row,:)=p{i}(i,:); p{i}(i,:)=mid;
和 不 同 的 右 端 向 量 b ， 给 定 迭 代 误 差 要 求 ， 用 Jacobi 和
Gauss-Seidel 迭代法计算，观测得出的迭代向量序列是否收敛。若收敛，记录迭代次数，分 析计算结果并得出你的结论。 (b) 选定初始向量初始向量 x
(0)
和不同的右端向量 b ，如取 x
(0)
= 0, b = Ae, e = (1,1,⋯1)T 。
n= 17 （2）取 x
(0)
= 0 ， b =（1,1， ……，1,1）T,迭代误差取默认值
10
−5
，迭代次数上限取默认值
50，使用 Gauss-Seidel 法进行迭代。 >> test2
>> b=ones(20,1) >> x0=zeros(20,1) >> [x,n]=Gauss_SeidelMethod(A,b,x0) x= 0.2766 0.2327 0.2159 0.2223 0.2228 0.2221 0.2222 0.2222 0.2222 0.2222 0.2222 0.2222 0.2222 0.2222 0.2221 0.2228 0.2223 0.2159 0.2327 0.2766
1.3 M 文件
.m 1.3.1 CreateMatrix CreateMatrix.m function [A,b]=CreateMatrix(n) %用于存放习题1的题目信息，并建立构造题目中矩阵的函数 %对矩阵A赋值 A1=6*ones(1,n); A2=ones(1,n-1); A3=8*ones(1,n-1); A=diag(A1)+diag(A2,1)+diag(A3,-1); %对向量b赋值 b=15*ones(n,1); b(1)=7; b(n)=14;
-0.0889 0.3111 -0.0889 0.3111 -0.0889 0.3111 -0.0889 0.3111 -0.0889 0.3111 -0.0882 0.3105 -0.0957 0.3313 -0.0472
n= 16 （4）取 x
(0)
= 0 ， b =（1,0,1,0,……，1,0）T,迭代误差取默认值
b= 7 15 15 15
15 15 15 15 15 14 >> x=GaussMethod(A,b) x= 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 n=84 时， >> [A,b]=CreateMatrix(84) >> x=GaussMethod(A,b) 1.0e+008 * 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 …… -0.0001 0.0002 -0.0003 0.0007 -0.0013 0.0026 -0.0052 0.0105
…… 1.0000 1.0000 1.0000 1.0000
1.2 计算结果分析
对于高斯消元法， 当 n=10 时， 得到的计算结果为 x= （1.0000， 1.0000， 1.0000， 1.0000， ……， T 1.0000） ， 恰好为精确值； 但当 n=84 时， 得到的计算结果却为 x=(0.0000,0.0000,0.0000,……， -0.0001,0.0002， ……，2.3487，-4.0263，5.3684）T,完全偏离了精确值。 对于列主元消去法，当 n=10 时，得到的计算结果为 x= （ 1.0000 ， 1.0000 ， 1.0000 ， T 1.0000， ……， 1.0000） ， 恰好为精确值； 当 n=10 时， 得到的计算结果为 x= （1.0000， 1.0000， 1.0000，1.0000， ……，1.0000）T，恰好为精确值。 由于使用高斯消元法时可能会出现小主元作除数的情况，虽然在求解系数矩阵为 10 阶 的方程时，由于步数较少，尚能求得精确解，但在求解求解系数矩阵为 84 阶的巨型方程时， 由于求解步数太多，使得小主元作除数的情况累加，最终导致“大数除小数”的情况发生， 因而所得解严重偏离精确解。 而使用列主元消去法时， 由于有效地避免了这一不利情况， 因而无论是在求解系数矩阵 为 10 阶还是 84 阶的方程时，都能得到精确解。