Chapter3-System of Equation

In this chapter, we consider the problem of solving several simultaneous Equations in several variables. Most of our attention will be paid to the case where the Number of equations and the number of unknown variables are the same.

Gaussian Elimination

We begin by describing the simplest form of Gaussian elimination. In fact, it is so Simple that it is not guaranteed to proceed to completion, let alone find an accurate Solution. The modifications that will be needed to improve the “naive” method will be Introduced beginning in the next section. Three useful operations can be applied to a linear system of equations that yield An equivalent system, meaning one that has the same solutions. These operations are As follows:

1. Swap one equation for another.
2. Add or subtract a multiple of one equation from another.
3. Multiply an equation by a nonzero constant.

Example of GE

Apply Gaussian elimination in tableau form for the system of three equations in three unknowns:

$$\begin{array}{c}x+2y-z=3\\ 2x+y-2z=3\\ -3x+y+z=-6.\end{array}$$

This is written in tableau form as

$$\left[\begin{array}{rrrr}1& 2& -1& 3\\ 2& 1& -2& 3\\ -3& 1& 1& -6\end{array}\right]$$

Two steps are needed to eliminate column 1: and one more step to eliminate column 2: Returning to the equations

$$\begin{array}{rl}x+2y-z& =3\\ -3y& =-3\\ -2z& =-4,\end{array}$$

we can solve for the variables

$$\begin{array}{rl}x& =3-2y+z\\ -3y& =-3\\ -2z& =-4\end{array}$$

and solve for $z,y,x$ in that order. The latter part is called back substitution, or back solving because, after elimination, the equations are readily solved from the bottom up. The solution is $x=3,y=1,z=2$ .

Generative Form

For system of equations like:

$$\left[\begin{array}{cccccc}{a}_{11}& a_{12}& \dots & a_{1n}& |& b_1\\ a_{21}& a_{22}& \dots & a_{2n}& |& b_2\\ ⋮& ⋮& & ⋮& |& ⋮\\ a_{n1}& a_{n2}& \dots & a_{nn}& |& b_n\end{array}\right]$$

We use these update-formula to update the matrix:

$$a_{ij}\leftarrow a_{ij}-\left(\frac{a_{i1}}{a_{11}}\right)a_{1j}{b}_i\leftarrow b_i-\left(\frac{a_{i1}}{a_{11}}\right)b_1,a_{11}\ne 0$$

After the elimination is completed, the tableau is upper triangular:

$$\left[\begin{array}{cccccc}{a}_{11}& a_{12}& \dots & a_{1n}& |& b_1\\ 0& a_{22}& \dots & a_{2n}& |& b_2\\ ⋮& ⋮& & ⋮& |& ⋮\\ 0& 0& \dots & a_{nn}& |& b_n\end{array}\right]$$

After elimination, the equations are readily solved from the bottom up using back solving.

$$\begin{array}{c}{x}_n=b_n/a_{nn}\\ x_{n-1}=\left(b_{n-1}-a_{n-1,n}{x}_n\right)/a_{n-1,n-1}\\ x_i=\left(b_i-\sum _{j=i+1}^n{a}_{ij}{x}_j\right)/a_{ii},i=n-1,n-2,\dots ,1\end{array}$$

LU Factorization

Carrying the idea of tableau form one step farther brings us to the matrix form of a system of equations. Matrix form will save time in the long run by simplifying the algorithms and their analysis.

The $LU$ factorization is a matrix representation of Gaussian elimination. It consists of writing the coefficient matrix $A$ as a product of a lower triangular matrix $L$ and an upper triangular matrix $U$. The $LU$ factorization is the Gaussian elimination version of a long tradition in science and engineering—breaking down a complicated object into simpler parts.

$$\left[\begin{array}{ccccc}{a}_{11}& a_{12}& a_{11}& ...& a_{11}\\ a_{21}& a_{22}& a_{22}& ...& a_{2n}\\ a_{31}& a_{32}& a_{32}& ...& a_{3n}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ a_{n2}& a_{n2}& a_{n2}& ...& a_{nn}\end{array}\right]=\left[\begin{array}{ccccc}{l}_{11}& 0& 0& ...& 0\\ l_{21}& l_{22}& 0& ...& 0\\ l_{31}& l_{32}& l_{53}& ...& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ l_{n1}& l_{n2}& l_{n3}& ...& l_{nn}\end{array}\right]\cdot \left[\begin{array}{ccccc}{u}_{11}& u_{12}& u_{13}& ...& u_{1n}\\ 0& u_{22}& u_{23}& ...& u_{2n}\\ 0& 0& u_{33}& ...& u_{2n}\\ 0& ⋮& ⋮& \ddots & u_{3n}\\ 0& 0& 0& ...& u_{nn}\end{array}\right]$$

• LU matrix An $m\times n$ matrix $L$ is lower triangular if its entries satisfy $l_{i,j}=0$ for $i<j$. An $m\times n$ matrix $U$ is upper triangular if its entries satisfy $u_{i,j}=0$ for $i>j$

Example of LU Factorization

Find the LU factorization of

$$A=\left[\begin{array}{rrr}1& 2& -1\\ 2& 1& -2\\ -3& 1& 1\end{array}\right]$$

This matrix is the matrix of coefficients of system (2.4). The elimination steps proceed as before:

$$\begin{array}{rl}& \left[\begin{array}{rrr}1& 2& -1\\ 2& 1& -2\\ -3& 1& 1\end{array}\right]⟶\begin{array}{r}\text{ subtract }2\times \text{ row }1\\ \text{ from row 2 }\end{array}⟶\left[\begin{array}{rrr}1& 2& -1\\ 0& -3& 0\\ -3& 1& 1\end{array}\right]\\ & ⟶\begin{array}{r}\text{ subtract }-3\times \text{ row }1\\ \text{ from row }3\end{array}⟶\left[\begin{array}{rrr}1& 2& -1\\ 0& -3& 0\\ 0& 7& -2\end{array}\right]\\ & ⟶\text{ from row }3⟶\left[\begin{array}{rrr}1& 2& -1\\ 0& -3& 0\\ 0& 0& -2\end{array}\right]=U\\ & \end{array}$$

The lower triangular $L$ matrix is formed, as in the previous example, by putting 1 ‘s on the main diagonal and the multipliers in the lower triangle-in the specific places they were used for elimination. That is,

$$L=\left[\begin{array}{rrr}1& 0& 0\\ 2& 1& 0\\ -3& -\frac{7}{3}& 1\end{array}\right]$$

Two Forms of LU Factorization

Doolittle Factorization

$$\left[\begin{array}{ccccc}{a}_{11}& a_{12}& a_{11}& ...& a_{11}\\ a_{21}& a_{22}& a_{22}& ...& a_{2n}\\ a_{31}& a_{32}& a_{32}& ...& a_{3n}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ a_{n2}& a_{n2}& a_{n2}& ...& a_{nn}\end{array}\right]=\left[\begin{array}{ccccc}1& 0& 0& ...& 0\\ l_{21}& 1& 0& ...& 0\\ l_{31}& l_{32}& 1& ...& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ l_{n1}& l_{n2}& l_{n3}& ...& 1\end{array}\right].\left[\begin{array}{ccccc}{u}_{11}& u_{12}& u_{13}& ...& u_{1n}\\ 0& u_{22}& u_{23}& ...& u_{2n}\\ 0& 0& u_{33}& ...& u_{2n}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& ...& u_{nn}\end{array}\right]$$

Groot Factorization

$$\left[\begin{array}{ccccc}{a}_{11}& a_{12}& a_{11}& ...& a_{11}\\ a_{21}& a_{22}& a_{22}& ...& a_{2n}\\ a_{31}& a_{32}& a_{32}& ...& a_{3n}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ a_{n2}& a_{n2}& a_{n2}& ...& a_{nn}\end{array}\right]=\left[\begin{array}{ccccc}{l}_{11}& 0& 0& \dots & 0\\ l_{21}& l_{22}& 0& \dots & 0\\ l_{31}& l_{32}& l_{33}& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ l_{n1}& l_{n2}& l_{n3}& \dots & l_{nn}\end{array}\right]\cdot \left[\begin{array}{ccccc}1& u_{12}& u_{13}& \dots & u_{1n}\\ 0& 1& u_{23}& \dots & u_{2n}\\ 0& 0& 1& \dots & u_{3n}\\ ⋮& ⋮& ⋮& \ddots & \\ 0& 0& 0& \dots & ⋮\end{array}\right]$$

The fact that not all matrices have an LU factorization means that more work is required before we can declare the LU factorization a general Gaussian elimination algorithm. The related problem of swamping is described in the next section. The $PA=LU$ factorization is introduced, which will overcome both problems.

Sources to Error

Let $x_a$ be an approximate solution of the linear system $Ax=b$ . The residual is the vector $r=b-A{x}_a$ . The backward error is the norm of the residual $\Vert b-A{x}_a{\Vert }_{\infty }$ , and the forward error is $\Vert x-x_a{\Vert }_{\infty }$ .( ${\Vert A\Vert }_{\infty }={\max }_{1\le i\le n}\left(\left|a_{i1}\right|+\cdots +\left|a_{in}\right|\right)$ )

Denote the residual by $r=b-A{x}_a$ . The relative backward error of system $Ax=b$ is defined to be

$$\frac{\Vert r{\Vert }_{\infty }}{\Vert b{\Vert }_{\infty }}=\frac{\Vert b-A{x}_a{\Vert }_{\infty }}{\Vert b{\Vert }_{\infty }}$$

and the relative forward error is

$$\frac{\Vert x-x_a{\Vert }_{\infty }}{\Vert x{\Vert }_{\infty }}$$

The error magnification factor for $Ax=b$ is the ratio of the two:

$$\begin{array}{rl}eooromagnificationfactor& =\frac{relativeforwarderror}{relativebackwarderror}\\ & =\frac{\frac{\Vert x-x_a{\Vert }_{\infty }}{\Vert x{\Vert }_{\infty }}}{\frac{\Vert b-A{x}_a{\Vert }_{\infty }}{\Vert b{\Vert }_{\infty }}}\end{array}$$

The condition number of a square matrix $A$ , $cond(A)$ , is the maximum possible error magnification factor for solving $Ax=b$ , over all right-hand sides $b$ . Namely:

$$cond\left(A\right)=\underset{b\in {\mathbb{R}}^n,b\ne 0}{\max }\frac{\Vert x_a-x{\Vert }_{\infty }/\Vert x{\Vert }_{\infty }}{\Vert r{\Vert }_{\infty }/\Vert b{\Vert }_{\infty }}$$

For an invertible matrix $A$ , its condition number is $cond(A)=\Vert A{\Vert }_{\infty }\cdot \Vert A^{-1}{\Vert }_{\infty }$

Swamping

A second significant source of error in classical Gaussian elimination is much easier to fix. We demonstrate swamping with the next example.

$$E.g.\{\begin{array}{c}10^{-2}{x}_1+x_2=1\\ x_1+2x_2=4\end{array},itssolution:x_2=\frac{2}{1-2\times 10^{-20}},x_2=\frac{1-4\times 10^{-20}}{1-2\times 10^{-2}}$$

Gaussian elimination should be kept as small as possible to avoid swamping. Fortunately, there is a simple modification to native Gaussian elimination that forces the absolute value of multipliers to be no larger than 1. This new protocol, which involves judicious row exchanges in the tableau, is called partial pivoting.

Partical pivoting gaussian elimination

Example of PPGE

Apply Gaussian elimination with partial pivoting to solve the system

$$\begin{array}{rl}{x}_1-x_2+3x_3& =-3\\ -x_1-2x_3& =1\\ 2x_1+2x_2+4x_3& =0.\end{array}$$

This example is written in tableau form as

$$\left[\begin{array}{rrrr}1& -1& 3& -3\\ -1& 0& -2& 1\\ 2& 2& 4& 0\end{array}\right]\text{. }$$

Under partial pivoting we compare $\left|a_{11}\right|=1$ with $\left|a_{21}\right|=1$ and $\left|a_{31}\right|=2$, and choose $a_{31}$ for the new pivot. This is achieved through an exchange of rows 1 and 3:

$$\begin{array}{rl}& \left[\begin{array}{rrrr}1& -1& 3& -3\\ -1& 0& -2& 1\\ 2& 2& 4& 0\end{array}\right]\to \begin{array}{c}\text{ exchange row 1 }\\ \text{ and row 3 }\end{array}⟶\left[\begin{array}{rrrr}2& 2& 4& 0\\ -1& 0& -2& 1\\ 1& -1& 3& -3\end{array}\right]\\ & ⟶\begin{array}{r}\text{ subtract }-\frac{1}{2}\times \text{ row }1\\ \text{ from row }2\end{array}⟶\left[\begin{array}{rrrr}2& 2& 4& 0\\ 0& 1& 0& 1\\ 1& -1& 3& -3\end{array}\right]\\ & ⟶\begin{array}{r}\text{ subtract }\frac{1}{2}\times \text{ row }1\\ \text{ from row }3\end{array}⟶\left[\begin{array}{rrrr}2& 2& 4& 0\\ 0& 1& 0& 1\\ 0& -2& 1& -3\end{array}\right]\\ & \end{array}$$

Before eliminating column 2 we must compare the current $\left|a_{22}\right|$ with the current $\left|a_{32}\right|$. Because the latter is larger, we again switch rows:

$$\begin{array}{rl}& \left[\begin{array}{rrrr}2& 2& 4& 0\\ 0& 1& 0& 1\\ 0& -2& 1& -3\end{array}\right]⟶\begin{array}{r}\text{ exchange row }2\\ \text{ and row }3\end{array}⟶\left[\begin{array}{rrrr}2& 2& 4& 0\\ 0& -2& 1& -3\\ 0& 1& 0& 1\end{array}\right]\\ & ⟶\text{ from row 3 }\to \left[\begin{array}{rrrr}2& 2& 4& 0\\ 0& -2& 1& -3\\ 0& 0& \frac{1}{2}& -\frac{1}{2}\end{array}\right]\\ & \end{array}$$

Note that all three multipliers are less than 1 in absolute value. The equations are now simple to solve. From

$$\begin{array}{rl}\frac{1}{2}{x}_3& =-\frac{1}{2}\\ -2x_2+x_3& =-3\\ 2x_1+2x_2+4x_3& =0,\end{array}$$

we find that $x=[1,1,-1]$.

PA = LU Factorization

we put together everything we know about Gaussian elimination into the PA = LU factorization. This is the matrix formulation of elimination with partial pivoting. The $PA=LU$ factorization is the established workhorse for solving systems of linear equations.

Example of find the PA=LU matrix of A

Example of PA=LU to solve system of equations

Use the $\mathrm{PA}=\mathrm{LU}$ factorization to solve the system $Ax=b$ , where

$$A=\left[\begin{array}{rrr}2& 1& 5\\ 4& 4& -4\\ 1& 3& 1\end{array}\right],b=\left[\begin{array}{l}5\\ 0\\ 6\end{array}\right].$$

It remains to complete the two back substitutions.

1. $Lc=Pb$ :

$$\left[\begin{array}{rrr}1& 0& 0\\ \frac{1}{4}& 1& 0\\ \frac{1}{2}& -\frac{1}{2}& 1\end{array}\right]\left[\begin{array}{l}{c}_1\\ c_2\\ c_3\end{array}\right]=\left[\begin{array}{lll}0& 1& 0\\ 0& 0& 1\\ 1& 0& 0\end{array}\right]\left[\begin{array}{l}5\\ 0\\ 6\end{array}\right]=\left[\begin{array}{l}0\\ 6\\ 5\end{array}\right]$$

Starting at the top, we have

$$\begin{array}{rl}{c}_1& =0\\ \frac{1}{4}(0)+c_2& =6\Rightarrow c_2=6\\ \frac{1}{2}(0)-\frac{1}{2}(6)+c_3& =5\Rightarrow c_3=8.\end{array}$$

2. $Ux=c:$

$$\left[\begin{array}{rrr}4& 4& -4\\ 0& 2& 2\\ 0& 0& 8\end{array}\right]\left[\begin{array}{l}{x}_1\\ x_2\\ x_3\end{array}\right]=\left[\begin{array}{l}0\\ 6\\ 8\end{array}\right]$$

Starting at the bottom,

$$\begin{array}{rl}8x_3=8& \Rightarrow x_3=1\\ 2x_2+2(1)=6& \Rightarrow x_2=2\\ 4x_1+4(2)-4(1)=0& \Rightarrow x_1=-1.\end{array}$$

Therefore, the solution is $x=[-1,2,1]$.

Iterative Methods

The generative form of iterative methods:

$$Q{x}_{k+1}=(Q-A)x_k+b$$

Where $Q$ is the Split Matrix. Rewrite it at fixed-point formula:

$$\begin{array}{rl}form:x_{k+1}& =(I-Q^{-1}A)x_k+Q^{-1}b\\ IterativeMatrix:G& =(I-Q^{-1}A)\end{array}$$

Richardson

$$Q{x}_{k+1}=(Q-A)x_k+b,Q=I,G=I-A$$

Example of Richardson iterative method.

$\{\begin{array}{c}3u+v=5\\ u+2v=5\end{array}$ , the coefficient matrix $A$ : $\left(\begin{array}{cc}3& 1\\ 1& 2\end{array}\right)$ Richardson iterative method:

$$\left[\begin{array}{c}{u}_{k+1}\\ v_{k+1}\end{array}\right]=\left[\begin{array}{cc}-2& -1\\ -1& -1\end{array}\right]\left[\begin{array}{c}{u}_k\\ v_k\end{array}\right]+\left[\begin{array}{c}5\\ 5\end{array}\right]$$

---

Jacobi

$$Q{x}_{k+1}=(Q-A)x_k+bQ=D,G=-D^{-1}(L+U)$$

Where $D$ is the diagonal matrix.

$$\begin{array}{c}D=\left[\begin{array}{cccc}{a}_{11}& 0& \cdots & 0\\ 0& a_{22}& \ddots & 0\\ ⋮& \ddots & \ddots & ⋮\\ 0& 0& \cdots & a_{nn}\end{array}\right]\\ \\ L=\left[\begin{array}{cccc}0& 0& \dots & 0\\ {\alpha }_{21}& 0& \dots & 0\\ ⋮& \ddots & \ddots & ⋮\\ {\alpha }_{n1}& {\alpha }_{n2}& \dots & 0\end{array}\right],U=\left[\begin{array}{cccc}0& a_{12}& \dots & a_{1n}\\ 0& 0& \dots & a_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & 0\end{array}\right]\end{array}$$

Example of Jacobi Iterative Method

$\{\begin{array}{c}3u+v=5\\ u+2v=5\end{array}$ , the correspond coefficient matrix $A$ : $\left(\begin{array}{cc}3& 1\\ 1& 2\end{array}\right)$ Jacobi Iterative method:

$$\begin{array}{rl}\left[\begin{array}{c}{u}_{k+1}\\ v_{k+1}\end{array}\right]& ={\left[\begin{array}{cc}3& 0\\ 0& 2\end{array}\right]}^{-1}\left[\begin{array}{cc}0& -1\\ -1& 0\end{array}\right]\left[\begin{array}{c}{u}_k\\ v_k\end{array}\right]+{\left[\begin{array}{cc}3& 0\\ 0& 2\end{array}\right]}^{-1}\left[\begin{array}{c}5\\ 5\end{array}\right]\\ & =\left[\begin{array}{cc}0& -1/3\\ -1/2& 0\end{array}\right]\left[\begin{array}{c}{u}_k\\ v_k\end{array}\right]+\left[\begin{array}{c}5/3\\ 5/2\end{array}\right]\\ & =\left[\begin{array}{c}(5-v_k)/3\\ (5-u_k)/2\end{array}\right]\end{array}$$

Relation between Jacobi and Richardson method

Gauss-Seidel

$$Q{x}_{k+1}=(Q-A)x_k+b,Q=D+L,G=-(L+D{)}^{-1}U$$

Example of Gauss-Seidel Method

$\{\begin{array}{c}3u+v=5\\ u+2v=5\end{array}$ , the coefficient matrix $A$ : $\left[\begin{array}{cc}3& 1\\ 1& 2\end{array}\right]$ Gauss-Seidel iterative method:

$$\left[\begin{array}{cc}3& 0\\ 1& 2\end{array}\right]\left[\begin{array}{c}{u}_{k+1}\\ v_{k+1}\end{array}\right]=\left[\begin{array}{cc}0& -1\\ 0& 0\end{array}\right]\left[\begin{array}{c}{u}_k\\ v_k\end{array}\right]+\left[\begin{array}{c}5\\ 5\end{array}\right]\Rightarrow \left[\begin{array}{c}{u}_{k+1}\\ v_{k+1}\end{array}\right]=\left[\begin{array}{c}\frac{5-v_k}{3}\\ \frac{5-u_{k+1}}{2}\end{array}\right]$$

Connection between Gauss-Seidel and Jacobi method

SOR

For generative linear iterative method:

$$x_{k+1}=G{x}_k+c$$

We can get a new iterative method by using relaxation method:

$$x_{k+1}=\gamma (G{x}_k+c)+(1-\gamma )x_k,\gamma \ne 0$$

Namely:

$$x_{k+1}=G_{\gamma }{x}_k+\gamma c,G_{\gamma }=\gamma G+(1-\gamma )I$$

SOR iterative method:

$$Q{x}_{k+1}=(Q-A)x_k+b,Q=\frac{1}{\omega }(D+\omega L)$$

Example of Gauss-Seidel Method

$\{\begin{array}{c}3u+v=5\\ u+2v=5\end{array}$ , the coefficient matrix $A$ : $\left[\begin{array}{cc}3& 1\\ 1& 2\end{array}\right]$ Gauss-Seidel iterative method:

SOR iterative method:

$$\begin{array}{rl}& \\ Gauss-Seidel:\left[\begin{array}{c}{u}_{k+1}\\ v_{k+1}\end{array}\right]& =\left[\begin{array}{c}\frac{5-v_k}{3}\\ \frac{5-u_{k+1}}{2}\end{array}\right]\\ SOR:\left[\begin{array}{c}{u}_{k+1}\\ v_{k+1}\end{array}\right]& =(1-\omega )\left[\begin{array}{c}{u}_k\\ v_k\end{array}\right]+\omega \left[\begin{array}{c}\frac{5-v_k}{3}\\ \frac{5-u_{k+1}}{2}\end{array}\right]\end{array}$$

Convergence of Iterative Methods

If the iterative matrix $G$ ‘s norm (of any kind) $\Vert G\Vert =\delta <1$ , then $x_{k+1}=(I-Q^{-1}A)x_k+Q^{-1}b$ is convergent to the solution of $Ax=b$ , also:

$$\Vert x-x_k\Vert \le \frac{\delta }{1-\delta }\Vert x_k-x_{k-1}\Vert \le \cdots \frac{{\delta }^k}{1-\delta }\Vert x_1-x_0\Vert$$

Above is exactly what strict diagonal Dominance implies

For $Ax=b$ , if $A$ is a strict diagonal matrix, then we have:

1. If $A$ is Dominance, Richardson iterative method is convergent.
2. If $A$ is strict row diagonal matrix, Jacobi iterative method is convergent.

Convergence of Richardson iterative method.

Convergence of Jacobi iterative method.

spectral radius method (spectral radius=$\sqrt{\max{\lambda_i}}$)

Convergence Judgment Using Spectral Radius method.

Example of Spectral Radius Method

fix: $\rho(G)>1$

Methods for Symmetric Positive-definite Matrices

• Definition(Symmetric positive-definite matrixes) The matrix $A\in {\mathbb{R}}^{n\times n}$ is symmetric if $A^{\top }=A$ ; A is positive-definite if $x^{\top }Ax>0$ for all vector $x\ne 0$.

If the $n\times n$ matrix $A$ is symmetric, then $A$ is positive-definite if and only if all of its eigenvalues are positive.

Property

• Symmetric positive-definite matrixes is nonsingular.
• If $A$ is symmetric positive-definite matrix, $P$ is a full-rank $n\times m$ matrix, then $P^{\top }AP$ is also a symmetric positive-definite matrix.

Cholesky Factorization

• Cholesky Factorization Theorem If $A$ is a symmetric positive-definite $n\times n$ matrix, then there exists an upper triangular $n\times n$ matrix $R$ such that $A=R^{\top }R$ .

Cholesky factorization

for k =1, 2, ... n
	if A(k,k) < 0, stop, end
	R(k,k)=math.sqrt(A(k,k))
	u=A(k, k+1:n)/R(k,k)
	R(k, k+1:n)=u
	A(k+1: n,k+1:n)-=np.dot(np.reverse(u),u)
end

Numerical Example

SOR Iterative Method of Symmetric Positive-definite Matrix

If $A$ is a symmetric positive-definite matrix, then the SOR iterative method for solving $Ax=b$ is convergent for any initial guess.

Conjugate Gradient Method

• A-Conjugate Let $A$ be a symmetric positive-definite $n\times n$ matrix. For two n-vectors $v$ and $w$, define the A-inner product $(v,w{)}_A=v^{\top }Aw$The vectors $v$ and $w$ are A-conjugate if $(v,w{)}_A=0$.

code of CGM