佛山科学技术学院
实 验 报 告
课程名称 数值分析
实验项目 数值积分
专业班级 机械工程 姓 名 余红杰 学 号 2111505010 指导教师 陈剑 成 绩 日 期 月 日
一、实验目的
1、 掌握程序的录入和matlab 的使用和操作；
2、 了解影响线性方程组解的精度的因素——方法与问题的性态。
3、 学会Matlab 提供的“\”的求解线性方程组。
二、实验要求
1、按照题目要求完成实验内容；
2、写出相应的Matlab 程序；
3、给出实验结果(可以用表格展示实验结果)；
4、分析和讨论实验结果并提出可能的优化实验。
5、写出实验报告。
三、实验步骤
1、用LU 分解及列主元高斯消去法解线性方程组
a)⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----15900001.582012151526099999.2310
7104321x x x x ， 输出b Ax =中系数LU A =分解的矩阵L 和U ，解向量x 和)det(A ;用列主元法的行交换次序解向量x 和求)det(A ；比较两种方法所得结果。
2、用列主高斯消元法解线性方程组b Ax =。
（1）、⎪⎪⎪⎭
⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎭⎫ ⎝⎛--11134.981.4987.023.116.427
.199.103.601.3321x x x
（2）、⎪⎪⎪⎭
⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎭⎫ ⎝⎛--11134.981.4990.023.116.427
.199.103.600.3321x x x 分别输出)det(,,A b A ，解向量x ，（1）中A 的条件数。
分析比较（1）、（2）的计算结果
3、线性方程组b Ax =的A 和b 分别为
⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=1095791068565778710A ,⎪⎪⎪⎪⎪⎭
⎫ ⎝⎛=31332332b 则解T x ),1,1,1,1(=. 用MATLAB 内部函数求)det(A 和A 的所有特征值和2)(A cond . 若令
⎪⎪⎪⎪⎪⎭
⎫ ⎝⎛=+98.99599.6989.998.585604.508.72.71.8710A A δ, 求解b x x A A =++))((δδ，输出向量x δ和2x δ，从理论结果和实际计算两方面分析线性
方程组b Ax =解的相对误差22/x x
δ以及A 的相对误差
/A A δ的关系。
四、实验结果
1：
%run311.m
clc,clear;
A = [10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];
b = [8;5.90001;5;1];
%L U 分解
format short
%小数点后四位，不然会受到后面的影响
[L U] = lu(A)
%解方程组，输出A ，det(A)
y = L\b;
format long
%小数点后15位显示
x = U\y
%V ol_xiao.m列主元消去，未用增广矩阵而把系数矩阵A和矩阵B分开对应来的function x=V ol_xiao(A,b)
%列主元消去
%x为解
%A为系数矩阵
n =length(A);
%A默认为方阵，求其大小，且默认和b长度一样
C =zeros(1,n);
%用于系数矩阵交换存放数据
b0 =0;
%用于b矩阵交换存放数据
x =zeros(n,1);
%建一个矩阵存放解
%下面对矩阵进行n-1次列主元交换得到上三角矩阵
for i=1:(n-1)
for j=(i+1):n
if abs(A(i,i))<abs(A(j,i))
C =A(i,1:n);
A(i,1:n)=A(j,1:n);
A(j,1:n)=C;
b0 = b(i);
b(i) =b(j);
b(j) =b0;
else continue;
end
end;
for j=(i+1):n
b(j) = b(j) + b(i)*(-A(j,i))/A(i,i);
A(j,1:n)=A(j,1:n)+A(i,1:n).*(-A(j,i)/A(i,i));
end
end
%以下为判断是否满秩然后从下往上求解：
if A(n,n)~=0
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
else x = ['err'];
end
%run312.m
clc,clear;
A = [10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];
b = [8;5.90001;5;1];
format long
x=V ol_xiao(A,b)
（其实可以看出，两种结果是一致的）
2.
%run321.m
clc;clear;
A = [3.01 6.03 1.99;1.27 4.16 -1.23;0.987 -4.81 9.34];
b = [1;1;1];
format bank
%格式转换为普通的，不然会用科学计数法表示
x = Vol_xiao(A,b)
T1 = cond(A,1)
T2 = cond(A,2)
T3 = cond(A,inf)
%可以看出条件数很大
%run322.m
clc;clear;
A = [3.01 6.03 1.99;1.27 4.16 -1.23;0.990 -4.81 9.34]
b = [1;1;1]
det1 = det(A)
format bank
x = Vol_xiao(A,b)
%从条件数就可以猜出其实微小变化就会对结果造成很大差异
%结果也验证了条件数大的方程的不稳定性很差。
3．
%run331.m
clc;clear;
A = [10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10];
b = [32;23;33;31];
format short
%显示对应行列式值和特征值
det1 = det(A)
eig1 = eig(A)
format bank
%显示条件数
tjs2 = cond(A,2)
%run332.m
clc;clear;
A = [10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10];
b = [32;23;33;31];
X = [1;1;1;1];
AA= [10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98];
oA= AA - A;
XX= AA \ b;
oX= XX – X
%输出oX 的2范数
normx2 = norm(oX)
%输出X 的相对误差
xrerro = normx2/norm(X)
%输出A 的相对误差
Arerro = norm(oA)/norm(A)
五、讨论分析
1.对于第一题，有一个顺序主子式是极小的，接近等于0，所以无法直接得到标准的上下三角，所以调用程序内置的lu 函数进行分解，所求到的值与列主元方法的值一致，查找后发现matlab 的lu 分解就是根据列主元方法求出的；
2.
条件数对一个方程的稳定性有很大影响，题中的系数矩阵条件数超过五万，微
小的变化就引起了结果的很大变化，如果是数据分析的话相当于蝴蝶效应般严重；
3.第三题也是对于一个系数矩阵的范数，很小的变化引起结果的大变化，该系数矩阵的相关数也是三千多，也是比较大，不稳定的。
六、改进实验建议
相关题目太少，而且很少有扩展性题目，对于各种指令的使用并不充分，其实在这个求解方程组的过程中，还是有比较多的相关方法可以进行MATLAB实现的。