数值分析（四）：解习题中的三两分乐趣

于我而言，做作业，解习题是件无聊，甚至是反感而近于厌恶的一件事情，尤其是面对一些作业量大，却又重复使用同一个理论进行证明的题目，既浪费精力，又搅乱心神，极有一种撕碎撒花，碾压脚踏载有习题的这几页书的冲动，但是这种近乎野蛮无赖的方式竟从没有发生在我身上，只能说在习题上燃起的聒噪的情绪，在闲书处又得到了安置和净化，每每有这种情绪，都选择避开这些魍魉，转向静谧的闲书中沉浸一会儿，静下心来，再回到教材中去升华自己的知识体系，以便铲除这些障碍，虽然最后肯定也是马马虎虎地交待了任务，但是情绪上的体验必然是不佳的，不过这次的习题中，倒有一点值得自己欣慰，因为完成习题的方法，毕竟不是小学中学那样，只能借助纸笔和迟钝的脑袋进行演算，而现在，习题的复杂性增加，相应的，解习题的工具也是丰富成熟了不少，道高一尺魔高一丈意义大抵如此，这个题目用到了C++，所以在这里记录下自己这段再也不会二次出现的感觉。
题目如下：

SOR法也不是什么高深的新方法，此刻忽然想起了炒概念这么一件事，很多人都以炒概念为能，吹得天花乱坠，以致让自己产生了一种能力强、且学识渊博的错觉，笔者曾经也中过毒，但是忘了是什么时候醒悟的，现在回想一通，还是非常汗颜，以后碰到炒概念的人，自己沉默即可，毕竟大家都不容易，人艰不拆。回到SOR法，这是一种迭代法，是高斯-赛德尔迭代法的一种改进，该方法认为高斯-赛德尔的那种使用最新更新的数据还是收敛得不够快，所以就想到用一个松弛因子，加速序列的收敛，于是就有了超松弛迭代法，超松弛迭代法的迭代公式如下：
既然是高斯-赛德尔迭代法的改进，所以第一步还是老样子：

然后松弛因子就派上用场了，将没有更新数据的xi(k)同上面这个高斯-赛德尔迭代法的结果做一个加权，使得加权之后的x能够更快收敛于精确解，然后作为迭代公式的下一次已知数据，公式如下：

为了更快收敛于精确解，往往这个加权之后的相对于xi(k)和高斯-赛德尔迭代算的那个x（k+1）是往外推的，也就是说松弛因子w一般是比1大一点的数字，这样才能更快收敛于精确解。
所以，根据SOR法，得到本题目的迭代公式如下：

题目要求是根据这么一个精度，这个是迭代法求得的近似解向量和精确解向量的差，这个差也是一个向量，精度要求是这个差向量的无穷范数小于5乘以十的负六次方，也就是这个差向量里面的元素的绝对值最大的那一个都要小于这个数字，就是这么一个度量。
题目呢就是这么一个题目，思路也是简单得不能再简单，但是就是有点畏难，你想万一要迭代个十几次，那我用计算器算岂不是真的极费心神，犹豫了一会，到底用不用计算器算，这就是当时把迭代公式写出来后的我的短暂的心路历程，后来一想还是不愿算。
既然是迭代，那就计算机来算吧，反正计算机就是干这行的，所以就用C++来帮我实现这一枯燥的过程了，解特定的题目，那就不要弄得太抽象，也不是针对一般性的线性方程，所以，干脆地就针对这个迭代公式，敲几行代码，很容易就把迭代次数给求出来了，接下来就是简单的思路：
解向量，直接用vector <double存储，和以往一样，采用函数的形式，输入和输出，输入是解向量（初始解向量就是0向量），松弛因子，以及误差要求，输出就是一个解向量，所以函数的形式如下：

vector<double> D(vector<double> a, double e, double w)

然后是函数内部，既然是迭代，按照迭代的规矩来，必然有迭代结束的条件，也就是误差的要求：

    vector<double> s(3);
	s[0] = a[0] - 0.5;
	s[1] = a[1] - 1;
	s[2] = a[2] + 0.5;
	double max = fabs(s[0]);
	for (double i : s)
		if (max < fabs(i))
			max = fabs(i);
	if (max < e)
		return a;

迭代向量和精确向量之间的差向量，这个差向量里面的元素绝对值最大的那个也要小于5乘以十的负六次方，这样才能满足要求；
然后是主要的迭代部分，既然迭代公式都有了，就简单明了l了，直接用代码，把这段迭代公式表示一下就OK了：

a[0] = (1 - w)*a[0] + (w / 4)*(a[1] + 1);
a[1] = (1 - w)*a[1] + (w / 4)*(a[0] + a[2] + 4);
a[2] = (1 - w)*a[2] + (w / 4)*(a[1] - 3);
n++;
cout << n << "  :";
for (double j : a)
	cout << j << " ";
	cout << endl;
	return D(a, e, w);

这样，专门针对这个迭代公式的递推就是这样了，知道满足了精度要求就会推出递推，这里的n是迭代次数，我把它设置了全局变量，接下来是整个代码：

#include<vector>
#include<iostream>
#include<cmath>

using namespace std;

int n = 0;


vector<double> D(vector<double> a, double e, double w)
{
	vector<double> s(3);
	s[0] = a[0] - 0.5;
	s[1] = a[1] - 1;
	s[2] = a[2] + 0.5;
	double max = fabs(s[0]);
	for (double i : s)
		if (max < fabs(i))
			max = fabs(i);
	if (max < e)
		return a;
	else
	{
	a[0] = (1 - w)*a[0] + (w / 4)*(a[1] + 1);
	a[1] = (1 - w)*a[1] + (w / 4)*(a[0] + a[2] + 4);
	a[2] = (1 - w)*a[2] + (w / 4)*(a[1] - 3);
	n++;
	cout << n << "  :";
	for (double j : a)
		cout << j << " ";
	cout << endl;
	return D(a, e, w);
	}
}

int main()
{
	double e = 0.000005, w1 = 1.03, w2 = 1, w3 = 1.1;
	vector<double> x{ 0,0,0 };
	cout << "当w=1.03时，迭代情况如下：" << endl;
	D(x, e, w1);

	cout << "当w=1时，迭代情况如下：" << endl;
	n = 0;
	D(x, e, w2);

	cout << "当w=1.1时，迭代情况如下：" << endl;
	n = 0;
	D(x, e, w3);

	system("pause");
	return 0;
}

然后是主函数调用结果：

所以，很清楚地看到每次递推的结果，明白地看到松弛因子分别取这些值的时候，要5,6,6次迭代才能达到相应的精度要求。
在写习题的时候，能够用C++帮着解决点问题，还是比较有意思的，尤其是最后我验证这个结果对不对的时候，发现真的是哭笑不得，为什么呢？因为自己手算的时候，代数都代错了，自己没有用到最新算的数据，结果沦落到了用的是雅克比迭代法，这样肯定是得不到正确答案的呀，所以，最后搞得自己有点呵呵了，虽然计算机只会做你让它去做的事情，但是你让它做的事情，只要命令没错，它就完成得比你快而且准确，这也算是做题当中的一点小乐趣吧，哈哈，继续加油！