🥥MATLAB与Python数值计算

使用直接方法求解线性系统

解线性系统方法

矩阵分解技术

迭代和最小二乘法解线性系统

迭代方法

雅可比法

给定方程的线性系统：

$\boldsymbol{a}{11}\boldsymbol{x}1+\boldsymbol{a}{12}\boldsymbol{x}2+\cdots +\boldsymbol{a}{1\boldsymbol{n}}\boldsymbol{x}{\boldsymbol{n}}=\boldsymbol{b}1\\\boldsymbol{a}{21}\boldsymbol{x}1+\boldsymbol{a}{22}\boldsymbol{x}2+\cdots +\boldsymbol{a}{2\boldsymbol{n}}\boldsymbol{x}{\boldsymbol{n}}=\boldsymbol{b}2\\\vdots \\\boldsymbol{a}{\boldsymbol{i}1}\boldsymbol{x}1+\boldsymbol{a}{\boldsymbol{i}2}\boldsymbol{x}2+\cdots +\boldsymbol{a}{\boldsymbol{in}}\boldsymbol{x}{\boldsymbol{n}}=\boldsymbol{b}{\boldsymbol{i}}\\\vdots \\\boldsymbol{a}{\boldsymbol{n}1}\boldsymbol{x}1+\boldsymbol{a}{\boldsymbol{n}2}\boldsymbol{x}2+\cdots +\boldsymbol{a}{\boldsymbol{nn}}\boldsymbol{x}{\boldsymbol{n}}=\boldsymbol{b}{\boldsymbol{b}}$

从上面的等式可以得出：

$\boldsymbol{x}1=\frac{1}{\boldsymbol{a}{11}}\left( \boldsymbol{b}1-\boldsymbol{a}{12}\boldsymbol{x}2-\cdots -\boldsymbol{a}{\boldsymbol{n}1}\boldsymbol{x}{\boldsymbol{n}} \right) \\\boldsymbol{x}2=\frac{1}{\boldsymbol{a}{22}}\left( \boldsymbol{b}2-\boldsymbol{a}{21}\boldsymbol{x}1-\boldsymbol{a}{23}\boldsymbol{x}3-\cdots -\boldsymbol{a}{\boldsymbol{n}1}\boldsymbol{x}{\boldsymbol{n}} \right) \\\vdots \,\, \vdots \\\boldsymbol{x}{\boldsymbol{i}}=\frac{1}{\boldsymbol{a}{\boldsymbol{ii}}}\left( \boldsymbol{b}{\boldsymbol{i}}-\boldsymbol{a}{\boldsymbol{i}1}\boldsymbol{x}1-\cdots -\boldsymbol{a}{\boldsymbol{ii}-1}\boldsymbol{x}{\boldsymbol{i}-1}-\boldsymbol{a}{\boldsymbol{ii}+1}\boldsymbol{x}{\boldsymbol{i}+1}-\cdots -\boldsymbol{a}{\boldsymbol{in}}\boldsymbol{x}{\boldsymbol{n}} \right) \\\vdots \,\, \vdots \\\boldsymbol{x}{\boldsymbol{n}}=\frac{1}{\boldsymbol{a}{\boldsymbol{nn}}}\left( \boldsymbol{b}{\boldsymbol{n}}-\boldsymbol{a}{\boldsymbol{n}1}\boldsymbol{x}1-\cdots -\boldsymbol{a}{\boldsymbol{nn}-1}\boldsymbol{x}{\boldsymbol{n}-1} \right)$

雅可比法是一种迭代方法，它从对解 $\left[ \boldsymbol{x}{1}^{\left( 0 \right)},\boldsymbol{x}{2}^{\left( 0 \right)},\cdots ,\boldsymbol{x}_{\boldsymbol{n}}^{\left( 0 \right)} \right] ^{\boldsymbol{T}}$ 的初始猜测开始。然后，使用迭代 $k$ 中的解来找到迭代 $k+1$ 中系统解的近似值。完成如下：

$\boldsymbol{x}{1}^{\left( \boldsymbol{k}+1 \right)}=\frac{1}{\boldsymbol{a}{11}}\left( \boldsymbol{b}1-\boldsymbol{a}{12}\boldsymbol{x}{2}^{\left( \boldsymbol{k} \right)}-\cdots -\boldsymbol{a}{\boldsymbol{n}1}\boldsymbol{x}{\boldsymbol{n}}^{\left( \boldsymbol{k} \right)} \right) \\\boldsymbol{x}{2}^{\left( \boldsymbol{k}+1 \right)}=\frac{1}{\boldsymbol{a}{22}}\left( \boldsymbol{b}2-\boldsymbol{a}{21}\boldsymbol{x}{1}^{\left( \boldsymbol{k} \right)}-\boldsymbol{a}{23}\boldsymbol{x}{3}^{\left( \boldsymbol{k} \right)}-\cdots -\boldsymbol{a}{\boldsymbol{n}1}\boldsymbol{x}{\boldsymbol{n}}^{\left( \boldsymbol{k} \right)} \right) \\\vdots \,\, \vdots \\\boldsymbol{x}{\boldsymbol{i}}^{\left( \boldsymbol{k}+1 \right)}=\frac{1}{\boldsymbol{a}{\boldsymbol{ii}}}\left( \boldsymbol{b}{\boldsymbol{i}}-\boldsymbol{a}{\boldsymbol{i}1}\boldsymbol{x}{1}^{\left( \boldsymbol{k} \right)}-\cdots -\boldsymbol{a}{\boldsymbol{ii}-1}\boldsymbol{x}{\boldsymbol{i}-1}^{\left( \boldsymbol{k} \right)}-\boldsymbol{a}{\boldsymbol{ii}+1}\boldsymbol{x}{\boldsymbol{i}+1}^{\left( \boldsymbol{k} \right)}-\cdots -\boldsymbol{a}{\boldsymbol{in}}\boldsymbol{x}{\boldsymbol{n}}^{\left( \boldsymbol{k} \right)} \right) \\\vdots \,\, \vdots \\\boldsymbol{x}{\boldsymbol{n}}^{\left( \boldsymbol{k}+1 \right)}=\frac{1}{\boldsymbol{a}{\boldsymbol{nn}}}\left( \boldsymbol{b}{\boldsymbol{n}}-\boldsymbol{a}{\boldsymbol{n}1}\boldsymbol{x}{1}^{\left( \boldsymbol{k} \right)}-\cdots -\boldsymbol{a}{\boldsymbol{nn}-1}\boldsymbol{x}{\boldsymbol{n}-1}^{\left( \boldsymbol{k} \right)} \right)$

通常，迭代 $\boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{k}+1 \right)}$ 中的解可以表示为：

$\boldsymbol{x}{\boldsymbol{i}}^{\left( \boldsymbol{k}+1 \right)}=\frac{1}{\boldsymbol{a}{\boldsymbol{ii}}}\left( \boldsymbol{b}{\boldsymbol{i}}-\sum{\boldsymbol{j}=1,\boldsymbol{j}\ne \boldsymbol{i}}^{\boldsymbol{n}}{\boldsymbol{a}{\boldsymbol{ij}}\boldsymbol{x}{\boldsymbol{j}}^{\left( \boldsymbol{k} \right)}} \right) ,\boldsymbol{i}=1,\cdots ,\boldsymbol{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)$

解决

我们写成：

$\boldsymbol{x}{1}^{\left( \boldsymbol{k}+1 \right)}=\frac{1}{2}\left( -1+\boldsymbol{x}{2}^{\left( \boldsymbol{k} \right)}-\boldsymbol{x}{3}^{\left( \boldsymbol{k} \right)} \right) \\\boldsymbol{x}{2}^{\left( \boldsymbol{k}+1 \right)}=\frac{1}{5}\left( 1+2\boldsymbol{x}{1}^{\left( \boldsymbol{k} \right)}+\boldsymbol{x}{3}^{\left( \boldsymbol{k} \right)} \right) \\\boldsymbol{x}{3}^{\left( \boldsymbol{k}+1 \right)}=\frac{1}{4}\left( 3-\boldsymbol{x}{1}^{\left( \boldsymbol{k} \right)}+2\boldsymbol{x}_{2}^{\left( \boldsymbol{k} \right)} \right)$

第一次迭代 $k=0$
$\boldsymbol{x}{1}^{\left( 1 \right)}=\frac{1}{2}\left( -1+\boldsymbol{x}{2}^{\left( 0 \right)}-\boldsymbol{x}{3}^{\left( 0 \right)} \right) =\frac{1}{2}\left( -1+0-0 \right) =-\frac{1}{2}\\\boldsymbol{x}{2}^{\left( 1 \right)}=\frac{1}{5}\left( 1+2\boldsymbol{x}{1}^{\left( 0 \right)}+\boldsymbol{x}{3}^{\left( 0 \right)} \right) =\frac{1}{5}\left( 1+2\left( 0 \right) +0 \right) =\frac{1}{5}\\\boldsymbol{x}{3}^{\left( 1 \right)}=\frac{1}{4}\left( 3-\boldsymbol{x}{1}^{\left( 0 \right)}+2\boldsymbol{x}_{2}^{\left( 0 \right)} \right) =\frac{1}{4}\left( 3-0+2\left( 0 \right) \right) =\frac{3}{4}$
第二次迭代 $k=1$
$\boldsymbol{x}{1}^{\left( 2 \right)}=\frac{1}{2}\left( -1+\boldsymbol{x}{2}^{\left( 1 \right)}-\boldsymbol{x}{3}^{\left( 1 \right)} \right) =\frac{1}{2}\left( -1+\frac{1}{5}-\frac{3}{4} \right) =-\frac{31}{40}\\\boldsymbol{x}{2}^{\left( 2 \right)}=\frac{1}{5}\left( 1+2\boldsymbol{x}{1}^{\left( 1 \right)}+\boldsymbol{x}{3}^{\left( 1 \right)} \right) =\frac{1}{5}\left( 1+2\cdot \frac{-1}{2}+\frac{3}{4} \right) =\frac{3}{20}\\\boldsymbol{x}{3}^{\left( 2 \right)}=\frac{1}{4}\left( 3-\boldsymbol{x}{1}^{\left( 1 \right)}+2\boldsymbol{x}_{2}^{\left( 1 \right)} \right) =\frac{1}{4}\left( 3-\frac{-1}{2}+2\frac{1}{5} \right) =\frac{39}{40}$
第三次迭代$k=2$
$\boldsymbol{x}{1}^{\left( 3 \right)}=\frac{1}{2}\left( -1+\boldsymbol{x}{2}^{\left( 2 \right)}-\boldsymbol{x}{3}^{\left( 2 \right)} \right) =\frac{1}{2}\left( -1+\frac{3}{20}-\frac{39}{40} \right) =-\frac{73}{80}\\\boldsymbol{x}{2}^{\left( 3 \right)}=\frac{1}{5}\left( 1+2\boldsymbol{x}{1}^{\left( 2 \right)}+\boldsymbol{x}{3}^{\left( 2 \right)} \right) =\frac{1}{5}\left( 1+2\cdot \frac{-31}{40}-\frac{-73}{80} \right) =\frac{17}{200}\\\boldsymbol{x}{3}^{\left( 3 \right)}=\frac{1}{4}\left( 3-\boldsymbol{x}{1}^{\left( 2 \right)}+2\boldsymbol{x}_{2}^{\left( 2 \right)} \right) =\frac{1}{4}\left( 3-\frac{31}{40}+2\frac{3}{20} \right) =\frac{163}{160}$

**示例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]]

最小二乘解

线性系统解中的病态和正则化技术

解非线性方程组

数据插补

微分和积分

解非线性常微分方程

常微分方程的非标准有限差分法

线性规划和二次规划

非线性规划

最优控制问题

🏈指点迷津 | Brief

PreviousPython模拟概率统计机器学习 NextPython对比MATLAB波动数据分析优势

Last updated 2 years ago

Was this helpful?