数值计算详细笔记（三）：非线性方程组解法

6. Linear Systems Ax = b

6.1 Basic Concepts

6.1.1 LSEs

the Linear System of Equations (LSEs):
$(I)\left\{ \begin{aligned} E_1: & a_{11}x_1+a_{12}x_2+...+a_{1n}x_n=b_1 \\ E_2: & a_{21}x_1+a_{22}x_2+...+a_{2n}x_n=b_2 \\ ... & \\ E_n: & a_{n1}x_1+a_{n2}x_2+...+a_{nn}x_n=b_n \end{aligned} \right.$(I)⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧​E1​:E2​:...En​:​a11​x1​+a12​x2​+...+a1n​xn​=b1​a21​x1​+a22​x2​+...+a2n​xn​=b2​an1​x1​+an2​x2​+...+ann​xn​=bn​​

6.1.2 Operations of LSEs

Multiplied operation - 数乘

Equation $E_i$Ei​ can be multiplied by any nonzero constant $\lambda$λ
$(\lambda E_i)\rightarrow E_i$(λEi​)→Ei​

Multiplied and added operation - 倍加

Equation $E_j$Ej​ can be multiplied by any nonzero constant $\lambda$λ, and added to Equation $E_i$Ei​ in place of $E_i$Ei​, denoted by
$(\lambda E_j+E_i)\rightarrow E_i$(λEj​+Ei​)→Ei​

Transposition - 交换

Equation $E_i$Ei​ and $E_j$Ej​ can be transposed in order, denoted by
$E_i \leftrightarrow E_j$Ei​↔Ej​

6.1.3 Augmented Matirx

$\tilde{A}=[A,\textbf{b}]= \left ( \begin{array}{c:c} \begin{matrix} a_{11}&a_{12}&...&a_{1n}\\ a_{21}&a_{22}&...&a_{2n}\\ ... & ... & ... &... \\ a_{n1}&a_{n2}&...&a_{nn}\\ \end{matrix}& \begin{matrix} b_1\\ b_2\\ ...\\ b_n \end{matrix} \end{array} \right )$A~=[A,b]=⎝⎜⎜⎛​a11​a21​...an1​​a12​a22​...an2​​............​a1n​a2n​...ann​​​b1​b2​...bn​​​⎠⎟⎟⎞​

6.2 Gaussian Elimination Method

6.2.1 Overall Description

The key point of Gaussian Elimination Method is changing the original matrix into upper-triangular matrix, then using backward–substitution method to calculate the answer.

6.2.2 Algorithm

• INPUT: $N$N-dimension, $A(N,N), B(N)$A(N,N),B(N)

• OUTPUT: Solution $x(N)$x(N) or Message that LESs has no unique solution.

• Step $1$1: For $k = 1,2,...,N-1$k=1,2,...,N−1, do step 2-4.

• Step $2$2: Set $p$p be the smallest integer with $k\leq p\leq N$k≤p≤N and $A_{p,k}\not= 0$Ap,k​​=0. If no $p$p can be found, output: “no unique solution exists”; stop.

• Step $3$3: If $p\not=k$p​=k, do transposition $E_p\leftrightarrow E_k$Ep​↔Ek​.

• Step $4$4: For $i=k+1,...,N$i=k+1,...,N

  1. Set $m_{i,k}=\displaystyle\frac{A(i,k)}{A(k,k)}$mi,k​=A(k,k)A(i,k)​
  2. Set $B(i)=B(i)-m_{i,k}B(k)$B(i)=B(i)−mi,k​B(k)
  3. For $j=k+1,...,N$j=k+1,...,N, set $A(i,j)=A(i,j)-m_{i,k}A(k,j)$A(i,j)=A(i,j)−mi,k​A(k,j);

• Step $5$5: If $A(N,N)\not=0$A(N,N)​=0, set $x(N)=\displaystyle\frac{B(N)}{A(N,N)}$x(N)=A(N,N)B(N)​; Else, output:“no unique solution exists.”

• Step $6$6: For $i=N-1,N-2,...,1$i=N−1,N−2,...,1, set
  $X(i)=[B(i)-\sum\limits_{j=i+1}^{N}A(i,j)x(j)]/A(i,i)$X(i)=[B(i)−j=i+1∑N​A(i,j)x(j)]/A(i,i)

• Step $7$7: Output the solution $x(N)$x(N).

6.3 Pivoting Strategies

6.3.1 Background

According to the process of Gaussian Elimination Method, We find that if $a_{kk}^{(k-1)}$akk(k−1)​ is too small, the roundoff error will be larger.
$m_{i,k}=\displaystyle\frac{A(i,k)}{A(k,k)}\\ X(i)=[B(i)-\sum\limits_{j=i+1}^{N}A(i,j)x(j)]/A(i,i)\\$mi,k​=A(k,k)A(i,k)​X(i)=[B(i)−j=i+1∑N​A(i,j)x(j)]/A(i,i)

Therefore, in order to reduce the roundoff error, we need to make the value of $a_{kk}^{(k-1)}$akk(k−1)​ larger.

6.3.2 Maximal Column Pivoting Technique

This method is to make $a_{kk}^{(k-1)}$akk(k−1)​ equal to the maximal value in its column.

6.3.3 Maximal Row Pivoting Technique

This method is to make $a_{kk}^{(k-1)}$akk(k−1)​ equal to the maximal value in its row.

6.3.4 Partial Pivoting Technique

This method is to make $a_{kk}^{(k-1)}$akk(k−1)​ equal to the maximal value in its remaining area.

6.3.5 Scaled Partial Pivoting Technique

$s_i=\max_{1\leq j\leq n}|a_{i,j}| \\ \displaystyle\frac{a_{kk}}{s_k}=\max_{k\leq i\leq n}\displaystyle\frac{a_{i,1}}{s_i}$si​=1≤j≤nmax​∣ai,j​∣sk​akk​​=k≤i≤nmax​si​ai,1​​

This method is to make $\displaystyle\frac{a_{kk}^{(k-1)}}{s_{k}}$sk​akk(k−1)​​ equal to the maximal value in its remaining area.

6.4 LU Factorization

6.4.1 The advantage of LU Factorization

$Ax=b\\ A=LU\\ L= \left( \begin{matrix} 1 & 0 & 0 & ... & 0 \\ l_{21} & 1 & 0 & ... & 0 \\ l_{31} & l_{32} & 1 & ... & 0 \\ ... & ... & ... & ... & ... \\ l_{n1} & l_{n2} & l_{n3} & ... & 1 \end{matrix} \right), R= \left( \begin{matrix} u_{11} & u_{12} & u_{13} & ... & u_{1n} \\ 0 & u_{22} & u_{23} & ... & u_{2n} \\ 0 & 0 & u_{33} & ... & u_{3n} \\ ... & ... & ... & ... & ... \\ 0 & 0 & 0 & ... & u_{nn} \end{matrix} \right)$Ax=bA=LUL=⎝⎜⎜⎜⎜⎛​1l21​l31​...ln1​​01l32​...ln2​​001...ln3​​...............​000...1​⎠⎟⎟⎟⎟⎞​,R=⎝⎜⎜⎜⎜⎛​u11​00...0​u12​u22​0...0​u13​u23​u33​...0​...............​u1n​u2n​u3n​...unn​​⎠⎟⎟⎟⎟⎞​

We can use two-step process to solve $LUx=b$LUx=b.

1. $y=Ux,Ly=b$y=Ux,Ly=b
2. Solve $Ly=b$Ly=b determining $y$y with forward substitution method.
3. Solve $Ux=y$Ux=y determining $x$x with forward substitution method.

6.4.2 LU Factorization through Gaussian Elimination

Theorem

If Gaussian elimination can be performed on the linear system $Ax=b$Ax=b without row interchanges, then the matrix $A$A can be factored into the product of a lower-triangular $L$L and an upper-triangular matrix $U$U,
$A=LU$A=LU
where
$L= \left( \begin{matrix} 1 & 0 & 0 & ... & 0 \\ m_{21} & 1 & 0 & ... & 0 \\ m_{31} & m_{32} & 1 & ... & 0 \\ ... & ... & ... & ... & ... \\ m_{n1} & m_{n2} & m_{n3} & ... & 1 \end{matrix} \right), R= \left( \begin{matrix} a_{11}^1 & a_{12}^1 & a_{13}^1 & ... & a_{1n}^1 \\ 0 & a_{22}^2 & a_{23}^2 & ... & a_{2n}^2 \\ 0 & 0 & a_{33}^3 & ... & a_{3n}^3 \\ ... & ... & ... & ... & ... \\ 0 & 0 & 0 & ... & a_{nn}^n \end{matrix} \right)$L=⎝⎜⎜⎜⎜⎛​1m21​m31​...mn1​​01m32​...mn2​​001...mn3​​...............​000...1​⎠⎟⎟⎟⎟⎞​,R=⎝⎜⎜⎜⎜⎛​a111​00...0​a121​a222​0...0​a131​a232​a333​...0​...............​a1n1​a2n2​a3n3​...annn​​⎠⎟⎟⎟⎟⎞​

Proof

$m_{j,1}=\displaystyle\frac{a_{j,1}}{a_{1,1}}\\ M^1 = \left( \begin{matrix} 1 & 0 & 0 & ... & 0 \\ -m_{21} & 1 & 0 & ... & 0 \\ -m_{31} & 0 & 1 & ... & 0 \\ ... & ... & ... & ... & ... \\ -m_{n1} & 0 & 0 & ... & 1 \end{matrix} \right)$mj,1​=a1,1​aj,1​​M1=⎝⎜⎜⎜⎜⎛​1−m21​−m31​...−mn1​​010...0​001...0​...............​000...1​⎠⎟⎟⎟⎟⎞​
Thus,
$A^n=M^{n-1}M^{n-2}...M^{1}A.$An=Mn−1Mn−2...M1A.
Let $U=A^n$U=An, then
$[M^1]^{-1}...[M^{n-2}]^{-1}[M^{n-1}]^{-1}U=A \\ [M^1]^{-1} = \left( \begin{matrix} 1 & 0 & 0 & ... & 0 \\ m_{21} & 1 & 0 & ... & 0 \\ m_{31} & 0 & 1 & ... & 0 \\ ... & ... & ... & ... & ... \\ m_{n1} & 0 & 0 & ... & 1 \end{matrix} \right)\\ L = [M^1]^{-1}...[M^{n-2}]^{-1}[M^{n-1}]^{-1}\\$[M1]−1...[Mn−2]−1[Mn−1]−1U=A[M1]−1=⎝⎜⎜⎜⎜⎛​1m21​m31​...mn1​​010...0​001...0​...............​000...1​⎠⎟⎟⎟⎟⎞​L=[M1]−1...[Mn−2]−1[Mn−1]−1

6.4.3 LU Factorization through Gaussian Elimination

$LU= \left( \begin{matrix} 1 & 0 & 0 & ... & 0 \\ l_{21} & 1 & 0 & ... & 0 \\ l_{31} & l_{32} & 1 & ... & 0 \\ ... & ... & ... & ... & ... \\ l_{n1} & l_{n2} & l_{n3} & ... & 1 \end{matrix} \right) \left( \begin{matrix} u_{11} & u_{12} & u_{13} & ... & u_{1n} \\ 0 & u_{22} & u_{23} & ... & u_{2n} \\ 0 & 0 & u_{33} & ... & u_{3n} \\ ... & ... & ... & ... & ... \\ 0 & 0 & 0 & ... & u_{nn} \end{matrix} \right)\\ LU=A= \left( \begin{matrix} a_{11} & a_{12} & a_{13} & ... & a_{1n} \\ a_{21} & a_{22} & a_{23} & ... & a_{2n} \\ a_{31} & a_{32} & a_{33} & ... & a_{3n} \\ ... & ... & ... & ... & ... \\ a_{n1} & a_{n2} & a_{n3} & ... & a_{nn} \end{matrix} \right)$LU=⎝⎜⎜⎜⎜⎛​1l21​l31​...ln1​​01l32​...ln2​​001...ln3​​...............​000...1​⎠⎟⎟⎟⎟⎞​⎝⎜⎜⎜⎜⎛​u11​00...0​u12​u22​0...0​u13​u23​u33​...0​...............​u1n​u2n​u3n​...unn​​⎠⎟⎟⎟⎟⎞​LU=A=⎝⎜⎜⎜⎜⎛​a11​a21​a31​...an1​​a12​a22​a32​...an2​​a13​a23​a33​...an3​​...............​a1n​a2n​a3n​...ann​​⎠⎟⎟⎟⎟⎞​

Algorithm

6.5 Strictly Diagonally dominant Matrix

6.5.1 Definition

The $n*n$n∗n matrix is said to be strictly diagonally dominant (严格对角占优) when
$|a_{ii}|>\sum\limits_{j=1,j\not=i}^{n} |a_{ij}|$∣aii​∣>j=1,j​=i∑n​∣aij​∣
holds for each $i=1,2,3,...,n$i=1,2,3,...,n.

6.5.2 Property

1. A strictly diagonally dominant matrix $A$A is nonsingular.
2. Moreover, in this case, Gaussian elimination can be performed on any linear system of the form $Ax=b$Ax=b to obtain its unique solution without row or column interchanges, and the computations are stable with respect to the growth of roundoff errors.

Proof for First Property

A matrix is singular means its determinant is zero.

A matrix’s determinant is zero means the n vectors in the matrix are linearly dependent.

Thus, matrix $A$A is singular means there exists a column vector $u$u that $Ax=0$Ax=0.

6.6 Positive Definite Symmetric Matrix

6.6.1 Definition

A matrix $A$A is positive definite if it is symmetric and if $x^TAx > 0$xTAx>0 for every $n$n-dimensional column vector $x\not=0$x​=0.

6.6.2 Property

If A is an $n*n$n∗n positive definite matrix, then

1. A is nonsingular;
2. $a_{ii}>0$aii​>0 for each $i=1,2,...,n$i=1,2,...,n;
3. $\max_{1\leq k,j\leq n}|a_{k,j}|\leq \max_{1\leq i\leq n}|a_{i,i}|$max1≤k,j≤n​∣ak,j​∣≤max1≤i≤n​∣ai,i​∣;
4. $a_{ij}^2<a_{ii}a_{jj}$aij2​<aii​ajj​ for each $i\not=j$i​=j.

6.6.3 Theorem

6.7 $LL^T$LLT Factorization

6.7.1 Definition

For a $n*n$n∗n symmetric and positive definite matrix $A$A with the form
$A= \left( \begin{matrix} a_{11} & a_{12} & a_{13} & ... & a_{1n} \\ a_{12} & a_{22} & a_{23} & ... & a_{2n} \\ a_{13} & a_{23} & a_{33} & ... & a_{3n} \\ ... & ... & ... & ... & ... \\ a_{1n} & a_{2n} & a_{3n} & ... & a_{nn} \end{matrix} \right)$A=⎝⎜⎜⎜⎜⎛​a11​a12​a13​...a1n​​a12​a22​a23​...a2n​​a13​a23​a33​...a3n​​...............​a1n​a2n​a3n​...ann​​⎠⎟⎟⎟⎟⎞​
where $A^T=A.$AT=A. We can factorize this matrix to the form like $LL^T=A$LLT=A, where $L$L is a lower triangular matrix with form as follows
$\left( \begin{matrix} l_{11} & 0 & 0 & \cdots & 0 \\ l_{21} & l_{22} & 0 & \cdots & 0 \\ l_{31} & l_{32} & l_{33} & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ l_{n1} & l_{n2} & l_{n3} & \cdots & l_{nn} \end{matrix} \right)$⎝⎜⎜⎜⎜⎜⎛​l11​l21​l31​⋮ln1​​0l22​l32​⋮ln2​​00l33​⋮ln3​​⋯⋯⋯⋱⋯​000⋮lnn​​⎠⎟⎟⎟⎟⎟⎞​

Thus, we need to determine the elements $l_{ij}$lij​, for $i\in[1,n]$i∈[1,n] and $j\in[1,n]$j∈[1,n].
$A= \left( \begin{matrix} a_{11} & a_{12} & a_{13} & ... & a_{1n} \\ a_{12} & a_{22} & a_{23} & ... & a_{2n} \\ a_{13} & a_{23} & a_{33} & ... & a_{3n} \\ ... & ... & ... & ... & ... \\ a_{1n} & a_{2n} & a_{3n} & ... & a_{nn} \end{matrix} \right)=LL^T$A=⎝⎜⎜⎜⎜⎛​a11​a12​a13​...a1n​​a12​a22​a23​...a2n​​a13​a23​a33​...a3n​​...............​a1n​a2n​a3n​...ann​​⎠⎟⎟⎟⎟⎞​=LLT

6.7.2 Choleski’s Algorithm

Calculate the value one row by one row

To factor the positive definite $n*n$n∗n matrix $A$A into $LL^T$LLT, where $L$L is lower triangular:

• INPUT: the dimension $n$n; entries $a_{ij}$aij​ of $A$A, for $i\in [1,n]$i∈[1,n] and $j\in[1,i]$j∈[1,i].

• OUTPUT: the entries $l_{ij}$lij​ of $L$L, for $i\in [1,n]$i∈[1,n] and $j\in[1,i]$j∈[1,i].

• Step $1$1: Set $l_{11} = \sqrt{a_{11}}$l11​=a11​​.

• Step $2$2: For $j\in[2,n]$j∈[2,n], set $l_{j1}=\displaystyle\frac{a_{1j}}{l_{11}}$lj1​=l11​a1j​​

• Step $3$3: For $i\in[2,n-1]$i∈[2,n−1], do Steps 4 and 5.

• Step $4$4: Set $l_{ii}=[a_{ii}-\sum_{j=1}^{i-1}l_{ij}^2]^{\frac{1}{2}}$lii​=[aii​−∑j=1i−1​lij2​]21​.

• Step $5$5: For $j\in[i+1,n]$j∈[i+1,n], set $l_{ji}=\displaystyle\frac{a_{ij}-\sum_{k=1}^{i-1}l_{ik}l_{jk}}{l_{ii}}$lji​=lii​aij​−∑k=1i−1​lik​ljk​​

• Step $6$6: Set $l_{nn}=[a_{nn}-\sum_{k=1}^{n-1}l_{nk}^2]^\frac{1}{2}$lnn​=[ann​−∑k=1n−1​lnk2​]21​

• Step $7$7: OUTPUT $l_{ij}$lij​ for $j\in[1,i]$j∈[1,i] and $i\in[1,n]$i∈[1,n]. STOP!

Example

6.8 $LDL^T$LDLT Factorization

6.8.1 Definition

Matrix $A$A is a positive definite matrix, thus
$A = LDL^T.$A=LDLT.

We can calculate the value of $L$L and $D$D one row by one row.

6.8.2 Algorithm

6.9 Tri-diagonal Linear System

6.9.1 Definition

An $n*n$n∗n matrix $A$A is called a band matrix (带状矩阵), if integers $p$p and $q$q, with $1<p, q<n$1<p,q<n, exist having the property that $a_{ij}=0$aij​=0 whenever $i+p\leq j$i+p≤j or $j+q\leq i$j+q≤i. The bandwidth (带宽) of a band matrix is defined as $w=p+q-1$w=p+q−1.

6.9.2 LU Factorization

$A = LU$A=LU

In order to solve the problem of $Ax=LUx=b$Ax=LUx=b, there are two steps to do.

1. $z = Ux$z=Ux, and solve $Lz=b$Lz=b
2. solve $Ux=z$Ux=z

6.9.3 Remarks

1. Band matrices usually are sparse matrices, thus we need to substitute two-dimensional array to one-dimensional array to store the value of the matrices.
2. Banded matrices appear in numerical calculation methods of partial differential equations and are common matrix forms.