使用直接方法求解线性系统 | 解线性系统方法 | 矩阵分解技术
迭代和最小二乘法解线性系统
迭代方法
雅可比法
给定方程的线性系统:
a11x1+a12x2+⋯+a1nxn=b1a21x1+a22x2+⋯+a2nxn=b2⋮ai1x1+ai2x2+⋯+ainxn=bi⋮an1x1+an2x2+⋯+annxn=bb
从上面的等式可以得出:
x1=a111(b1−a12x2−⋯−an1xn)x2=a221(b2−a21x1−a23x3−⋯−an1xn)⋮⋮xi=aii1(bi−ai1x1−⋯−aii−1xi−1−aii+1xi+1−⋯−ainxn)⋮⋮xn=ann1(bn−an1x1−⋯−ann−1xn−1)
雅可比法是一种迭代方法,它从对解[x1(0),x2(0),⋯,xn(0)]T的初始猜测开始。 然后,使用迭代k中的解来找到迭代k+1中系统解的近似值。 完成如下:
x1(k+1)=a111(b1−a12x2(k)−⋯−an1xn(k))x2(k+1)=a221(b2−a21x1(k)−a23x3(k)−⋯−an1xn(k))⋮⋮xi(k+1)=aii1(bi−ai1x1(k)−⋯−aii−1xi−1(k)−aii+1xi+1(k)−⋯−ainxn(k))⋮⋮xn(k+1)=ann1(bn−an1x1(k)−⋯−ann−1xn−1(k))
通常,迭代xi(k+1)中的解可以表示为:
xi(k+1)=aii1(bi−∑j=1,j=inaijxj(k)),i=1,⋯,n
对于任意\boldsymbol{\varepsilon }>0$,当$\left\| \boldsymbol{x}{\boldsymbol{i}}^{\left( \boldsymbol{k}+1 \right)}-\boldsymbol{x}{\boldsymbol{i}}^{\left( \boldsymbol{k} \right)} \right\| <\boldsymbol{\varepsilon }时,雅可比迭代停止。
**示例1,**为线性系统编写雅可比方法的前三个迭代:
$$ \left( \begin{matrix} 2& -1& 1\\ -2& 5& -1\\ 1& -2& 4\\\end{matrix} \right) \left( \begin{array}{c} \boldsymbol{x}_1\\ \boldsymbol{x}_2\\ \boldsymbol{x}_3\\\end{array} \right) =\left( \begin{array}{c} -1\\ 1\\ 3\\\end{array} \right) $$
从零向量开始
$$ \boldsymbol{x}^{\left( 0 \right)}=\left( \begin{array}{c} 0\\ 0\\ 0\\\end{array} \right) $$
解决
我们写成:
x1(k+1)=21(−1+x2(k)−x3(k))x2(k+1)=51(1+2x1(k)+x3(k))x3(k+1)=41(3−x1(k)+2x2(k))
第一次迭代k=0
x1(1)=21(−1+x2(0)−x3(0))=21(−1+0−0)=−21x2(1)=51(1+2x1(0)+x3(0))=51(1+2(0)+0)=51x3(1)=41(3−x1(0)+2x2(0))=41(3−0+2(0))=43
第二次迭代k=1
x1(2)=21(−1+x2(1)−x3(1))=21(−1+51−43)=−4031x2(2)=51(1+2x1(1)+x3(1))=51(1+2⋅2−1+43)=203x3(2)=41(3−x1(1)+2x2(1))=41(3−2−1+251)=4039
第三次迭代k=2
x1(3)=21(−1+x2(2)−x3(2))=21(−1+203−4039)=−8073x2(3)=51(1+2x1(2)+x3(2))=51(1+2⋅40−31−80−73)=20017x3(3)=41(3−x1(2)+2x2(2))=41(3−4031+2203)=160163
**示例2,**雅可比方法将用于求解线性系统:
$$ \left( \begin{matrix} -5& 1& -2\\ 1& 6& 3\\ 2& -1& -4\\\end{matrix} \right) \left( \begin{array}{c} \boldsymbol{x}_1\\ \boldsymbol{x}_2\\ \boldsymbol{x}_3\\\end{array} \right) =\left( \begin{array}{c} 13\\ 1\\ -1\\\end{array} \right) $$
MATLAB代码
function x = JacobiSolve(A, b, Eps)
n = length(b) ;
x0 = zeros(3, 1) ;
x = ones(size(x0)) ;
while norm(x-x0, inf) ≥ Eps
x0 = x ;
for i = 1 : n
x(i) = b(i) ;
for j = 1 : n
if j , i
x(i) = x(i) - A(i, j)*x0(j) ;
end
end
x(i) = x(i) / A(i, i) ;
end
>> A = [-5 1 -2; 1 6 3; 2 -1 -4] ;
>> b = [13; 1; -1] ;
>> JacobiSolveLinSystem(A, b)
-2.0000
1.0000
-1.0000
使用Python,功能JacobiSolve具有以下代码:
def JacobiSolve(A, b, Eps):
import numpy as np
n = len(b)
x0, x = np.zeros((n, 1), 'float'), np.ones((n, 1), 'float')
while np.linalg.norm(x-x0, np.inf) ≥ Eps:
x0 = x.copy()
for i in range(n):
x[i] = b[i]
for j in range(n):
if j != i:
x[i] -= A[i][j]*x0[j]
x[i] /= A[i][i]
return x
通过调用函数JacobiSolve求解给定的线性系统,我们获得:
In [3]: A = np.array([[-5, 1, -2], [1, 6, 3], [2, -1, -4]])
In [4]: b = np.array([[13], [1], [-1]])
In [5]: x = JacobiSolve(A, b, Eps)
In [6]: print(’x = \\n’, x)
Out[6]:
x =
[[-2. ],
[ 1.00000002],
[-1.00000002]]
最小二乘解
线性系统解中的病态和正则化技术
解非线性方程组
数据插补
微分和积分
解非线性常微分方程
常微分方程的非标准有限差分法
线性规划和二次规划
非线性规划
最优控制问题