中国矿业大学（北京）理学院
数值线性代数实验报告实验名称
组
长
一、实验目的，内容
四、数值结果平方根法求解方程组组员二、相关背景知识介绍五、计算结果的分析三、代码六、计算中出现的问题，解决方法及体会实验时间2014年10月9日
一、实验目的，内容
了解平方根法的原理和主要思想；
掌握平方根法的计算步骤；
用平方根法求解方程组，并观察其精度；
观察Hilbert 方程组的特点。
二、相关背景知识介绍
背景
解n 阶线性方程组Ax=b的choleskly 方法也叫做平方根法，这里对系数矩阵A 是有要求的，需要A 是对称正定矩阵，根据数值分析的相关理论，如果A 对称正定，那么系数矩阵就可以被分解为的A=L L T 形式，其中L 是下三角矩阵，将其代入Ax=b中，可得：LL T x=b进行如下分解：
y L T x L y b
那么就可先计算y, 再计算x ，由于L 是下三角矩阵，是L T 上三角矩阵，这样的计算比直接使用A 计算简便，同时你应该也发现了工作量就转移到了矩阵的分解上面，那么对于对称正定矩阵A 进行Cholesky 分解，其过程为：
设A=L L T ，即
a 11a 12 a 21a 22
a n 1a n 2... a 1n l 11 l 11l 21... l n 1 l ... a 2n l l ... l 22n 2 2122 ... ... ... ... ... ... ... a nn l n 1l n 2... l nn l nn
其中a ij a ji , i , j 1, 2,..., n
三、代码
问题：用平方根法解下列方程组，对计算解和准确解比较,观察准确程度
1 1
2（1） 1 3
1
4 1 [***********][1**********]181 5 1 x 1 b 1 6 x 2 b 2 1 ，右端项自己构造。 x 3 b 3 7 1 x 4 b 4 8 x 5 b 5 1 9
1 0. 500
（2） 0. 333 0. 250 0. 200 0. 5000. 3330. 2500. 200 x 1 b 1 0. 3330. 2500. 2000. 167 x 2 b 2 0. 2500. 2000. 1670. 143 x 3 b 3 ，再与（1）中的结果进行比较，你们 0. 2000. 1670. 1430. 125 x 4 b 4 0. 1670. 1430. 1250. 111 x 5 b 5
观察到什么现象？
[**************]先令x (1, 1, 1, 1, 1) T ，带入（1）中求得b (, , , , ) ，再将b 做为右端项[1**********]520
求解方程组。
a=[1,1/2,1/3,1/4,1/5;1/2,1/3,1/4,1/5,1/6;
1/3,1/4,1/5,1/6,1/7;1/4,1/5,1/6,1/7,1/8;1/5,1/6,1/7,1/8,1/9];%%输入a 矩阵
b=[137/6087/60153/140743/8401879/2520];%%输入b 矩阵
n=5;
l(1,1)=sqrt(a(1,1));
for i=2:1:n
l(i,1)=a(i,1)/l(1,1);
end %%计算矩阵l 的第一列元素for k=2:1:n
sum=0;
for p=1:1:k-1
sum=sum+l(k,p)*l(k,p);
end
l(k,k)=sqrt(a(k,k)-sum);%%计算矩阵l 的对角线元素
for i=k+1:1:n
sums=0;
for p=1:1:k-1
sums=sums+l(i,p)*l(k,p);
end
l(i,k)=(a(i,k)-sums)/l(k,k);
end
end %%计算矩阵l 的其他元素b(1)=b(1)/l(1,1);
for i=2:1:n
sum=0;
for j=1:1:i-1
sum=sum+l(i,j)*b(j);
end
b(i)=(b(i)-sum)/l(i,i);
end
b(n)=b(n)/l(n,n);
for i=n-1:-1:1
sum=0;
for j=i+1:1:n
sum=sum+l(j,i)*b(j);
end
b(i)=(b(i)-sum)/l(i,i);
end
b %%回代过程
四、数值结果
精确解
计算解
扰动后的计算解x (1, 1, 1, 1, 1) T x (1.0000,1.0000,1.0000,1.0000,1. 0000) T x (0.9784,1.2107,0.5170,1.3745,0. 9192) T
五、计算结果的分析
由计算结果可明显看出：（1）中的结果非常准确，但（2）的结果却与准确值之间的误差很大。然而，当方程组中，若当A 或者b 的只发生微小的变化而引起方程组解的很大变化，则称此方程组为病态方程组，A 称为病态矩阵。本题的方程组正好就是病态方程组。
用平方根法求解正定矩阵线性方程组的优点有以下两点：1）数值稳定。2）计算量小，大约为次乘除法，是一般矩阵A 的LU 分解计算量的一半。这样会导致误差增大且计算不易。
改进的平方根法就是引进中间量使得，这样就省去了开方运算且数值稳定精度大大的提高了。
六、计算中出现的问题，解决方法及体会
问题：
(1)编程计算矩阵A 的Cholesky 分解矩阵L 的几个for 循环出现紊乱，导致程序无法运行。
(2)求得矩阵L 之后进行最后一部计算求未知数x 时分别对i 和j 进行循环，j 循环中l(j,i)错写成l(i,j)从而导致无法运行出正确的结果。
解决方法：
(1)经过仔细推敲书本的知识，并且和身边的同学进行讨论得出正确的循环结构从而运行出正确的结果。
(2)在草稿纸进行演算，明白其深层含义从而纠正程序运行得到正确的结果。
体会：
(1)当for 循环结构较复杂时，应该回归到课本仔细推敲其算法的含义，逐步进行分析思考出正确的结果。
(2)上课的时候老师给出算法之后不仅仅只是照抄在书本上，应该仔细推敲其含义，不能因为抄错了算法从而一直无法运行出正确的结果，应该多思考正确答案。
教
师
评
语
指导教师：日期：年月日
中国矿业大学（北京）理学院
数值线性代数实验报告实验名称
组
长
一、实验目的，内容
四、数值结果平方根法求解方程组组员二、相关背景知识介绍五、计算结果的分析三、代码六、计算中出现的问题，解决方法及体会实验时间2014年10月9日
一、实验目的，内容
了解平方根法的原理和主要思想；
掌握平方根法的计算步骤；
用平方根法求解方程组，并观察其精度；
观察Hilbert 方程组的特点。
二、相关背景知识介绍
背景
解n 阶线性方程组Ax=b的choleskly 方法也叫做平方根法，这里对系数矩阵A 是有要求的，需要A 是对称正定矩阵，根据数值分析的相关理论，如果A 对称正定，那么系数矩阵就可以被分解为的A=L L T 形式，其中L 是下三角矩阵，将其代入Ax=b中，可得：LL T x=b进行如下分解：
y L T x L y b
那么就可先计算y, 再计算x ，由于L 是下三角矩阵，是L T 上三角矩阵，这样的计算比直接使用A 计算简便，同时你应该也发现了工作量就转移到了矩阵的分解上面，那么对于对称正定矩阵A 进行Cholesky 分解，其过程为：
设A=L L T ，即
a 11a 12 a 21a 22
a n 1a n 2... a 1n l 11 l 11l 21... l n 1 l ... a 2n l l ... l 22n 2 2122 ... ... ... ... ... ... ... a nn l n 1l n 2... l nn l nn
其中a ij a ji , i , j 1, 2,..., n
三、代码
问题：用平方根法解下列方程组，对计算解和准确解比较,观察准确程度
1 1
2（1） 1 3
1
4 1 [***********][1**********]181 5 1 x 1 b 1 6 x 2 b 2 1 ，右端项自己构造。 x 3 b 3 7 1 x 4 b 4 8 x 5 b 5 1 9
1 0. 500
（2） 0. 333 0. 250 0. 200 0. 5000. 3330. 2500. 200 x 1 b 1 0. 3330. 2500. 2000. 167 x 2 b 2 0. 2500. 2000. 1670. 143 x 3 b 3 ，再与（1）中的结果进行比较，你们 0. 2000. 1670. 1430. 125 x 4 b 4 0. 1670. 1430. 1250. 111 x 5 b 5
观察到什么现象？
[**************]先令x (1, 1, 1, 1, 1) T ，带入（1）中求得b (, , , , ) ，再将b 做为右端项[1**********]520
求解方程组。
a=[1,1/2,1/3,1/4,1/5;1/2,1/3,1/4,1/5,1/6;
1/3,1/4,1/5,1/6,1/7;1/4,1/5,1/6,1/7,1/8;1/5,1/6,1/7,1/8,1/9];%%输入a 矩阵
b=[137/6087/60153/140743/8401879/2520];%%输入b 矩阵
n=5;
l(1,1)=sqrt(a(1,1));
for i=2:1:n
l(i,1)=a(i,1)/l(1,1);
end %%计算矩阵l 的第一列元素for k=2:1:n
sum=0;
for p=1:1:k-1
sum=sum+l(k,p)*l(k,p);
end
l(k,k)=sqrt(a(k,k)-sum);%%计算矩阵l 的对角线元素
for i=k+1:1:n
sums=0;
for p=1:1:k-1
sums=sums+l(i,p)*l(k,p);
end
l(i,k)=(a(i,k)-sums)/l(k,k);
end
end %%计算矩阵l 的其他元素b(1)=b(1)/l(1,1);
for i=2:1:n
sum=0;
for j=1:1:i-1
sum=sum+l(i,j)*b(j);
end
b(i)=(b(i)-sum)/l(i,i);
end
b(n)=b(n)/l(n,n);
for i=n-1:-1:1
sum=0;
for j=i+1:1:n
sum=sum+l(j,i)*b(j);
end
b(i)=(b(i)-sum)/l(i,i);
end
b %%回代过程
四、数值结果
精确解
计算解
扰动后的计算解x (1, 1, 1, 1, 1) T x (1.0000,1.0000,1.0000,1.0000,1. 0000) T x (0.9784,1.2107,0.5170,1.3745,0. 9192) T
五、计算结果的分析
由计算结果可明显看出：（1）中的结果非常准确，但（2）的结果却与准确值之间的误差很大。然而，当方程组中，若当A 或者b 的只发生微小的变化而引起方程组解的很大变化，则称此方程组为病态方程组，A 称为病态矩阵。本题的方程组正好就是病态方程组。
用平方根法求解正定矩阵线性方程组的优点有以下两点：1）数值稳定。2）计算量小，大约为次乘除法，是一般矩阵A 的LU 分解计算量的一半。这样会导致误差增大且计算不易。
改进的平方根法就是引进中间量使得，这样就省去了开方运算且数值稳定精度大大的提高了。
六、计算中出现的问题，解决方法及体会
问题：
(1)编程计算矩阵A 的Cholesky 分解矩阵L 的几个for 循环出现紊乱，导致程序无法运行。
(2)求得矩阵L 之后进行最后一部计算求未知数x 时分别对i 和j 进行循环，j 循环中l(j,i)错写成l(i,j)从而导致无法运行出正确的结果。
解决方法：
(1)经过仔细推敲书本的知识，并且和身边的同学进行讨论得出正确的循环结构从而运行出正确的结果。
(2)在草稿纸进行演算，明白其深层含义从而纠正程序运行得到正确的结果。
体会：
(1)当for 循环结构较复杂时，应该回归到课本仔细推敲其算法的含义，逐步进行分析思考出正确的结果。
(2)上课的时候老师给出算法之后不仅仅只是照抄在书本上，应该仔细推敲其含义，不能因为抄错了算法从而一直无法运行出正确的结果，应该多思考正确答案。
教
师
评
语
指导教师：日期：年月日