---
title: 數學軟體實作 - Poisson's equation
tags: 2020 Fall - 數學軟體實作
GA: G-77TT93X4N1
slideOptions:
theme: beige
progress: true
controls: true
slideNumber: false
---

<style>
.reveal {
font-size: 24px;
}
</style>

# Lecture - Jacobi method to solve Poisson's equation

### References:
* [Poisson's equation](https://en.wikipedia.org/wiki/Poisson%27s_equation)
* [Jacobi method](https://en.wikipedia.org/wiki/Jacobi_method)

---

## Jacobi method - Recall

#### 解 $Ax=b$:

將原問題改寫成:
$$
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$.

---

所以 Jacobi method 最基本的樣子就是
<font color="red">
$$
x^{(n+1)}= D^{-1}\left(-(L+U)x^{(n)}+b\right)
$$
</font> 另外我們也會改寫上式為
$$
x^{(n+1)}= D^{-1}\left(Dx^{(n)}+b-(D+L+U)x^{(n)}\right) = x_n + D^{-1}(b-Ax^{(n)})
$$
接著我們定義 <font color="blue">$r^{(n)} = b-Ax^{(n)}$</font>, 這就是所謂的 residual, 那我們就有以下這個式子
<font color="blue">
$$
x^{(n+1)} = x^{(n)}+D^{-1}r^{(n)}
$$
</font> 依此我們可以定義兩種 stopping criteria:
* $\|x^{(n+1)}-x^{(n)}\|<Tol_1$
* $\|r^{(n)}\|<Tol_2$

---

## Poisson's equation in 1D

一維的 Poisson's equation 是以下形式:
$$
u''(x) = f(x), \quad x\in[0, 1],
$$
其中 $f(x)$ 是已知函數, 而 $u=u(x)$ 則是未知的函數, 通常需要伴隨一組邊界條件來確保解的唯一性, 例如 zero Dirichlet 邊界條件: $u(0)=u(1)=0$.

#### Example

例如, $f(x)=\pi^2\sin(\pi x)$, 也就是說想要找一個函數 $u(x)$ 使得它滿足
$$
u''(x) = \pi^2\sin(\pi x), \quad x\in[0, 1],
$$
並且 $u(0)=u(1)=0$.

好消息! 答案是 $u(x) = -\sin(\pi x)$. :100:

---

### Finite difference method

對於一般的 $f(x)$ 也許我們無法推導將問題解出來, 那就可以試著將原問題離散化求數值解.

意思就是, 原本我們是要求 $u(x)$, $x\in[0,1]$, 不過現在我們先將 $[0,1]$ 分成很多點
$$
x_k = k\Delta x, \quad k=1, \cdots, N, \quad \Delta x = \frac{1}{N+1}.
$$
並將 $u(x)$ 在這些點上的值定為 $u_k = u(x_k)$, $k=1,\cdots,N$.

 

---

然後我們試著將 $\{u_k\}$ 所要滿足的方程式寫下來. 而要滿足的方程式為以下這個:
$$
\frac{u_{k-1} - 2u_k + u_{k+1}}{(\Delta x)^2} = f_k, \quad k=1, \cdots, N, \quad u_0=u_{n+1}=0
$$

或寫成矩陣的形式:

$$
\frac{1}{(\Delta x)^2}
\begin{bmatrix}
-2 & 1 & & & & \\
1 & -2 & 1 & & & \\
& 1 & \ddots & \ddots & & \\
& & \ddots & \ddots & 1 & \\
& & & 1 & -2
\end{bmatrix}_{N\times N} u = f,
$$
其中
$$
u = [u_1, u_2, \cdots, u_N]^T, \quad f = [f_1, f_2, \cdots, f_N]^T.
$$
也就是說我們把 "解一維 Poisson's equation" 這件事轉換成"解一個線性系統".

---

接著我們就可以用 Jacobi method 來解它.

Given $x^{(n)}$,
$$
r^{(n)} = b-Ax^{(n)}, \quad x^{(n+1)} = x^{(n)}+D^{-1}r^{(n)}
$$

---

我們先看一下 Poisson 問題的 residual 怎麼算:

$$
\frac{1}{(\Delta x)^2}
\begin{bmatrix}
-2 & 1 & & & & \\
1 & -2 & 1 & & & \\
& 1 & \ddots & \ddots & & \\
& & \ddots & \ddots & 1 & \\
& & & 1 & -2
\end{bmatrix}_{N\times N} u = f,
$$
給定 $u^{(n)}$, residual 就是
$$
r^{(n)}_k = f_k - \frac{1}{(\Delta x)^2}\left(u^{(n)}_{k-1} - 2u^{(n)}_k + u^{(n)}_{k+1}\right)
$$
然後下一個 $u^{(n+1)}$ 就是 ($u^{(n+1)} = u^{(n)}+D^{-1}r^{(n)}$)
$$
u^{(n+1)} = u^{(n)}-\frac{(\Delta x)^2}{2}r^{(n)}
$$

---

### Summary

Given: $f(x)$

1. Compute $x_k$ and $f_k$
2. Initialize $u^{(0)}_k=0$
3. For $n\ge 1$
1. $r^{(n)}_k = f_k - \frac{1}{(\Delta x)^2}\left(u^{(n)}_{k-1} - 2u^{(n)}_k + u^{(n)}_{k+1}\right)$
2. $u^{(n+1)} = u^{(n)}-\frac{(\Delta x)^2}{2}r^{(n)}$

stopping criteria:
* $\|u^{(n+1)}-u^{(n)}\|<Tol_1$
* $\|r^{(n)}\|<Tol_2$

---

### Assignment

以 Jacobi method 求解 1D Poisson's equation
$$
u''(x) = f(x), \quad x\in[0, 1],
$$
邊界條件: $u(0)=u(1)=0$ 並且 $f(x)=\pi^2\sin(\pi x)$.

---

## Poisson's equation in 2D

二維的 Poisson's equation 是以下形式:
$$
u_{xx}+u_{yy} = f(x,y), \quad x\in[0, 1], \quad y\in[0,1]
$$
其中 $f(x,y)$ 是已知函數, 而 $u=u(x,y)$ 則是未知的函數, 通常需要伴隨一組邊界條件來確保解的唯一性, 例如 zero Dirichlet 邊界條件.

#### Example

例如, $f(x,y)=2\pi^2\sin(\pi x)\sin(\pi y)$, 也就是說想要找一個函數 $u(x,y)$ 使得它滿足
$$
u_{xx}+u_{yy} = 2\pi^2\sin(\pi x)\sin(\pi y), \quad x\in[0, 1], \quad y\in[0,1]
$$
並且 $u$ 在邊界為 $0$.

好消息! 答案是 $u(x,y) = -\sin(\pi x)*\sin(\pi y)$. :100:

---

### Finite difference method

二維的網格需要多加一個 $y$, 假設我們在 $x$ 方面切出 $N$ 個點, $y$ 方向切出 $M$ 個點:
$$
x_k = k\Delta x, \quad k=1, \cdots, N, \quad \Delta x = \frac{1}{N+1}, \\
y_l = l\Delta y, \quad l=1, \cdots, M, \quad \Delta y = \frac{1}{M+1}, \\
$$
並將 $u(x,y)$ 在這些點上的值定為 $u_{k,l} = u(x_k, y_l)$.


---

然後我們試著將 $\{u_{k,l}\}$ 所要滿足的方程式寫下來. 而要滿足的方程式為以下這個:
$$
\frac{u_{k-1,l} - 2u_{k,l} + u_{k+1,l}}{(\Delta x)^2}+\frac{u_{k,l-1} - 2u_{k,l} + u_{k,l+1}}{(\Delta y)^2} = f_{k,l}
$$

或寫成矩陣的形式: (for example) (注意以下正負號跟我們相反)
 or 

這是個 $NM\times NM$ 的線性系統

---

我們用 Jacobi method 來解它.
我們先看一下 Poisson 問題的 residual 怎麼算:

$$
\frac{u_{k-1,l} - 2u_{k,l} + u_{k+1,l}}{(\Delta x)^2}+\frac{u_{k,l-1} - 2u_{k,l} + u_{k,l+1}}{(\Delta y)^2} = f_{k,l}
$$
給定 $u^{(n)}$, residual 就是
$$
r^{(n)}_{k,l} = f_{k,l} - \left(\frac{u_{k-1,l} - 2u_{k,l} + u_{k+1,l}}{(\Delta x)^2}+\frac{u_{k,l-1} - 2u_{k,l} + u_{k,l+1}}{(\Delta y)^2}\right)
$$
然後下一個 $u^{(n+1)}$ 就是 ($u^{(n+1)} = u^{(n)}+D^{-1}r^{(n)}$)
$$
u^{(n+1)} = u^{(n)}-\frac{1}{2}\frac{1}{\frac{1}{(\Delta x)^2}+\frac{1}{(\Delta y)^2}}r^{(n)}
$$
如果 $\Delta x=\Delta y$, 則有 $u^{(n+1)} = u^{(n)}-\frac{(\Delta x)^2}{4}r^{(n)}$

---

### Summary

Given: $f(x,y)$, assuming that $N=M$ so $\Delta x=\Delta y$

1. Compute $x_k$, $y_l$ and $f_{k,l}$
2. Initialize $u^{(0)}_{k,l}=0$
3. For $n\ge 1$
1. $r^{(n)}_{k,l} = f_{k,l} - \left(\frac{u_{k-1,l} - 2u_{k,l} + u_{k+1,l}}{(\Delta x)^2}+\frac{u_{k,l-1} - 2u_{k,l} + u_{k,l+1}}{(\Delta y)^2}\right)$
2. $u^{(n+1)} = u^{(n)}-\frac{(\Delta x)^2}{4}r^{(n)}$

stopping criteria:
* $\|u^{(n+1)}-u^{(n)}\|<Tol_1$
* $\|r^{(n)}\|<Tol_2$

---

### Assignment

以 Jacobi method 求解 2D Poisson's equation
$$
u_{xx}+u_{yy} = f(x,y), \quad x\in[0, 1], \quad y\in[0,1]
$$
$u$ 在邊界都是 $0$ 並且 $f(x,y)=2\pi^2\sin(\pi x)\sin(\pi y)$.

---