线性代数方程组的数值解法

数学实验

高斯消元法

求解方程组$Ax=b$

• c++高斯消元模板:待补充
• matlab命令:反斜杠\
  
  • 若$A$为可逆方阵,输出原方程的解$x$,若$A$为$n\times m$矩阵$(n>m)$,这是一个超定方程组,且$A^TA$可逆,输出原方程的最小二乘解$x$.
• mathematica命令

A = {{1, 2}, {3, 4}}
B = {5, 6}
LinearSolve[A, B]
>> {-4, 9/2}

LU分解

维基百科

LU分解在本质上是高斯消元法的一种表达形式。实质上是将A通过初等行变换变成一个上三角矩阵，其变换矩阵就是一个单位下三角矩阵。这正是所谓的杜尔里特算法（Doolittle algorithm）：从下至上地对矩阵A做初等行变换，将对角线左下方的元素变成零，然后再证明这些行变换的效果等同于左乘一系列单位下三角矩阵，这一系列单位下三角矩阵的乘积的逆就是L矩阵，它也是一个单位下三角矩阵。

应用:

求解方程组$Ax=b$

1. 将A进行LU分解,$Ax=LUx=b$
2. 求解方程组$Ly=b$得到$y$
3. 求解方程组$Ux=y$得到$x$

注记:在两次的求解中，我们遇到的都是三角矩阵，因此运用向前（向后）替代法就可以简洁地求解（参见三角矩阵），而不需要用到高斯消元法。然而，在将A进行LU分解时，仍然要用到高斯消元法。因此，这个方法适合在要对许多个不同的b求解时用。

解的误差分析

由于实际问题导出的线性方程组$Ax=b$中,$A与b$往往带有误差(扰动),解得误差分析是讨论$A或b$的微小变化对$x$的影响.如果对系数的微小扰动引起解的很大变化,我们称为病态线性方程组,线性矩阵称为病态矩阵.
首先引入矩阵的范数
维基百科:矩阵范数
通过阅读wiki,我们发现只要是矩阵到实数的一个映射,满足严格正定性,线性性和三角不等式,就可称作一个范数
矩阵范数的度量有很多种,比较常见的是:

• p-范数诱导的矩阵范数
  

  当p = 2（欧几里德範数）且m = n（方陣）時，诱导的矩阵範数就是谱範数。矩陣A的谱範数是A最大的奇異值或半正定矩阵A*A的最大特徵值的平方根

求解线性方程组迭代方法(p100)

(p111 题5)已知方程组$Ax=b$,其中$A\in R^{20*20}$,是一个五对角矩阵,主对角线值均为$n$(自己定),旁边两条副对角线值为$-\frac{1}{2}$,再旁边两条值为$-\frac{1}{4}$,接下来我们采用两种方法,尝试不同的初值和$n$来迭代求解.

• 矩阵构造

n = 3;
B = diag(1/6*ones(1,19),1)+diag(1/6*ones(1,19),-1)+diag(1/12*ones(1,18),2)+diag(1/12*ones(1,18),-2)
%diag(a,i)代表主对角线向上偏移i条对角线为a向量,其他为0所生成的矩阵
%坑爹啊,matlab不能随便换行

• 雅克比迭代
  迭代公式为$x^{(k+1)}=B*x^{(k)}+b$,初值向量$x^{(0)}$自己选取,$B$矩阵可以由上面的$A$矩阵得到.

t = 10              %迭代次数
b = zero(20,1)      %20个零的列向量
x = rand(20,t);     %生成随机数作为初值
for i=2:t
    x(:,i)=A1*x(:,i-1)+b;
end
x                   %观察结果,发现快速收敛到0,正确
%我们采用以下两种不同的方程组右端项
b = ones(20,1)      %观察结果,发现快速收敛
b = [9/4,7/4,3/2*ones(1,16),7/4,9/4]
                    %注意此时b要转置再除以3
                    %这是凑好了解向量全部为1的情况,观察结果,发现快速收敛
%我们尝试修改初值的范围,并计算多少次迭代可以使迭代误差小于1e-5
x = 100*rand(20,t)
x = 100000*rand(20,t)
for i = 2:t
    if max(abs(x(:,i)-x(:,i-1)))<1e-5
       i
       break
    end
end
计算得33次迭代可以小于1e-5(指数速度收敛)
最后我们研究主对角线元素大小的影响
n = 10000;
format long
A1 = diag(1/(n*2)*ones(1,19),1)+diag(1/(n*2)*ones(1,19),-1)+diag(1/(n*4)*ones(1,18),2)+diag(1/(n*4)*ones(1,18),-2)
t = 50;
b = [9/4,7/4,3/2*ones(1,16),7/4,9/4]+n-3
x = 1*rand(20,t);
for i=2:t
    x(:,i)=A1*x(:,i-1)+b'/n;
end
x
for i = 2:t
    if max(abs(x(:,i)-x(:,i-1)))<1e-5
       i
       break
    end
end
发现n越大,收敛速度越快10-8,100-6,1000000-4

• 高斯-赛德尔迭代
  与雅克比迭代的区别仅仅在于解向量在迭代的过程中实时更新,比如计算$x^{(k+1)}_3$时利用$x^{(k+1)}_2$而不是$x^{(k)}_2$.
  以下是代码实现,分析过程同上

n = 3;
format long
A1= diag(1/6*ones(1,19),1)+diag(1/6*ones(1,19),-1)+diag(1/12*ones(1,18),2)+diag(1/12*ones(1,18),-2)
t = 50;
b = zeros(20,1);
x = 10000*rand(20,1);
for i=2:t
    for j=1:20
        x(j)=A1(j,:)*x+b(j);
    end
end
x