---
title: 數學軟體實作 - Solving linear system
tags: 2020 Fall - 數學軟體實作
GA: G-77TT93X4N1
---

# Lecture - Jacobi method to solve linear system

## 解線性聯立方程

### Reading:
* [張智星 - MATLAB程式設計：進階篇 -6-4 聯立線性方程式](http://mirlab.org/jang/books/matlabProgramming4guru/06-4_linear_algebra.asp?title=6-4%20%C1p%A5%DF%BDu%A9%CA%A4%E8%B5{%A6%A1)

### Rules
1. 解線性聯立方程 $Ax=b$, 使用 `A\b`
2. 求某反矩陣乘以某向量 $A^{-1}b$, 使用 `A\b`
3. **絕對禁止**求反矩陣
* **絕對禁止**`inv(A)*b` 這樣的做法
* 即使數學公式看起來需要求反矩陣, 也**絕對禁止**求反矩陣
* 例如將某向量 $b$ 投影到某平面, 公式為 $A(A^TA)^{-1}A^Tb$, 也**絕對不可以**這樣寫

            p = A*inv(A^T*A)*A^T*b

正確寫法是

            x = (A^T*A)\(A^T*b)
            p = A*x

---

## sequence of vectors

之前我們討論過數列, $\{a_n\}$, 不過這裡的 $a_n\in R$ 也就是一個數字. 現在我們要推廣一點看看如果 $a_n\in R^N$ 的時候要怎樣來看這數列, 具體來說就是我們想知道這數列會不會"收斂".

> 要談收斂之前需要先談"距離", 我們這邊直接討論更高級一點的 "norm".

### Vector norm

假設 $v = (v_1, v_2, \cdots, v_N)\in R^N$, 我們可以定義以下三種 norm:
$$
\| v\|_1 = |v_1| + |v_2| + \cdots + |v_N|
$$
$$
\| v\|_2 = \sqrt{v^2_1+v^2_2 + \cdots v^2_N}
$$
$$
\| v\|_{\infty} = \max_{1\le i\le N} |v_i|
$$
在 `matlab` 裡面, 一個向量的 norm 可以用 `norm` 這個指令來計算, 分別是

norm(V,2),    norm(V,1),   norm(V,Inf)

此外, 更廣義一點我們可以定義 p-norm:
$$
\| v\|_p = \left(|v_1|^p+|v_2|^p + \cdots |v_N|^p \right)^{1/p}
$$
這也可以用 `norm` 這個指令來算.

### Exercise

我想要知道 $\| x\|_p = 1$ 這個集合在二維長什麼樣子, 也就是以下這個集合:
$$
\{x \quad |\quad x\in R^2, \quad \| x\|_p = 1\}
$$

#### 想法:

我們知道
1. $\| x\|_2 = 1$ 是個單位圓, 所以 $p=1$ 的圖形是知道的
2. 單位向量就是長度等於 $1$ 的向量, 而任何一個向量只要把它除以自己的長度, 就變成單位向量了.

#### 做法:

* Input: 一個正數 p
* Output: p-norm 等於一的圖
* Step 1: 隨機找 N 個二維的點
* Step 2: 把每個點投影到 p-norm=1 的圖上
* Step 3: 畫出這些點

#### pseudo code

Input: 一個正數 p
Output: p-norm 等於一的圖
x= 隨機找 N 個二維的點
for ii=1:N
    x(ii) = x(ii)/norm(x(ii), p);
end
plot(x, '.')

---

### Limit of vector sequence

有了 norm 我們就可量距離, 那就把之前一模一樣的想法套用進來:

**已知 $\{a_n\}$ 是個收斂數列**, 欲判斷數列極限
$$
\lim_{n\to\infty} a_n = L.
$$
想法就是看 $n$ 很大的時候 $a_n$ 的值是多少, 如果沒有變動太多我們就覺得他收斂了!

也就是說, 我們會設定一個容許值(tolerance), 如果 $\|a_{n+1}-a_{n}\|$ 小於這個數字, 我們就覺得差不多了, 並且覺得 $L \approx a_{n+1}$.

---

### [Jacobi method](https://en.wikipedia.org/wiki/Jacobi_method)

給定一個方陣 $A$ 以及一個向量 $b$, 我們想要解出一個向量 $x$ 滿足 $Ax=b$.

#### 想法

將矩陣 $A$ 改寫成 $A=D+L+U$, 其中 $D$ 是對角線矩陣, 它保留 $A$ 的對角線部分並將其他非對角線部分都變成 $0$. 而 $L+U$ 則是把 $A$ 對角線部分變 $0$ 其他不變.

因此我們就可將原問題改寫成:
$$
Ax = b \quad\to\quad
(D+L+U)x=b
\quad\to\quad
Dx= -(L+U)x+b.
$$
然後我們可以定義一個數列
$$
Dx_{n+1}= -(L+U)x_{n}+b,
$$
或是
$$
x_{n+1}= D^{-1}\left(-(L+U)x_{n}+b\right),
$$
並且知道如果 $x_n$ 這個數列收斂了, 令他收斂到 $x$, 那 $x$ 就會滿足 $Ax=b$.

#### 做法 (example):

* Input: square matrix A and a vector b
* Output: a vector x satisfying Ax=b
* Step 1: split A into A = D + R, where D is the diagonal part (a vector), R is the rest
* Step 2: choose a random initial condition x_0
* Step 3: iterate until converge
* step 3.1: Evaluate y = -R*x_0 +b
* step 3.2: Solve Dx = y directly

#### pseudo code

Input: square matrix A and a vector b
Output: a vector x satisfying Ax=b

Step 1: obtain the diagonal vector d and the rest (also store as A)
for ii=1:N
    d(ii) = A(ii,ii)
    A(ii,ii)=0
end

Step 2: random initial condition
x0 = rand(N,1);
x1 = x0:

ite = 0;
while dif > Tol
    y = -A*x0 + b
    for ii=1:N
        x1(ii) = y(ii)/d(ii)
    end
    
    dif = distance between x0 and x1
    x0 = x1
    
    ite = ite+1
    if(ite > ite_max)
        warning('too many iteration')
        break
    end
    
end

x = x1