遇到隐式格式，我们需要求解一个线性方程组。怎么办呢？当然是Thomas algorithm

---

注意到第一处箭头。
接着是矩阵化：

对角占有就可以LU分解

---

重头戏：Thomas algorithm

$$A=\begin{bmatrix} b_1 &-c_1 & \\ -a_2 &b_2 &-c_2 \\ & \ddots\ & \ddots &\ddots\ \\ & & -a_{N-1} &b_{N-1}&c_{N-1}\\ & & &-a_{N} &b_{N} \end{bmatrix},L=\begin{bmatrix} 1 & & \\ e_2 &1 & \\ & \ddots\ & \ddots \\ & & e_{N-1} &1\\ & & &e_N &1 \end{bmatrix}$$
$$U=\begin{bmatrix} f &-c_1 & \\ &f_2 &-c_2 \\ & & \ddots &\ddots \\ & & &f_{N-1}&-c_{N-1}\\ & & &&f_N \end{bmatrix}$$
1.首先是计算 $f_1$
L的第一行乘以 U的第一列，得到 A的第一行第一列的元素，即 $$1*f_1=b_1$$
2.接着计算出每个 $e_i,f_i,i=2,3...N$
计算 $e_2$:
L的第二行乘以 U的第一列，得到 A第二行第一列的元素，即 $$e_2*f_1=-a_2$$
得到 $e_2=-a_2/f_1$
计算 $e_3$:
L的第三行乘以 U的第二列（乘以第一列为零），得到 A第三行第二列的元素，即 $$0*-c_1+e_3*f_2=-a_3$$
得到 $e_3=-a_3/f_2$
依次类推~~~
3.计算 $f_i,i=2,3..N$
L的第二行乘以 U的第二列，得到 A的第二行第二列的元素，即 $$e_2*-c_1+f_2=b_2$$
得到 $f_2=b_2+e_2*c_1$
计算 $f_3$:
L的第三行乘以 U的第三列，得到 A的第三行第三列的元素，即 $$e_3*-c_2+f_3=b_3$$
得到 $f_3=b_3+e_3*c_2$

依次类推

---

计算$y_i,i=1,2,3,...,N$

$$U=\begin{bmatrix} f &-c_1 & \\ &f_2 &-c_2 \\ & & \ddots &\ddots \\ & & &f_{N-1}&-c_{N-1}\\ & & &&f_N \end{bmatrix},X=\begin{bmatrix} x_1\\ x_2\\ \vdots\\ \vdots\\ x_N\\ \end{bmatrix}$$

$$Y=\begin{bmatrix} f &-c_1 & \\ &f_2 &-c_2 \\ & & \ddots &\ddots \\ & & &f_{N-1}&-c_{N-1}\\ & & &&f_N \end{bmatrix}*\begin{bmatrix} x_1\\ x_2\\ \vdots\\ \vdots\\ x_N\\ \end{bmatrix}=\begin{bmatrix} f_1*x_1-c_1x_2\\ f_2*x_2-c_2x_3\\ \vdots\\ \vdots\\ f_N*x_N\\ \end{bmatrix}=\begin{bmatrix} y_1\\ y_2\\ \vdots\\ \vdots\\ y_N\\ \end{bmatrix}$$

$$d=\begin{bmatrix} d_1\\ d_2\\ \vdots\\ \vdots\\ d_N\\ \end{bmatrix}$$

由 L的 $$L=\begin{bmatrix} 1 & & \\ e_2 &1 & \\ & \ddots\ & \ddots \\ & & e_{N-1} &1\\ & & &e_N &1 \end{bmatrix}$$
$L*Y=d$
L的第一行乘以Y，即 $1*y_1=d_1$, $y_1=d_1$
L的第二行乘以Y，即 $e_2*y_1+1*y_2=d_2$, $y_2=d_2-e_2*y_1$
L的第三行乘以Y，即 $e_3*y_2+1*y_3=d_3$, $y_3=d_2-e_2*y_2$
依次类推，
$y_i = d_i -d_i*y_{i-1},i=2,3,...,N$

---

计算$x_i,i=N,N-1,...,1$
由上面计算出$y_i$

$$Y=\begin{bmatrix} f &-c_1 & \\ &f_2 &-c_2 \\ & & \ddots &\ddots \\ & & &f_{N-1}&-c_{N-1}\\ & & &&f_N \end{bmatrix}*\begin{bmatrix} x_1\\ x_2\\ \vdots\\ \vdots\\ x_N\\ \end{bmatrix}=\begin{bmatrix} f_1*x_1-c_1x_2\\ f_2*x_2-c_2x_3\\ \vdots\\ \vdots\\ f_N*x_N\\ \end{bmatrix}=\begin{bmatrix} y_1\\ y_2\\ \vdots\\ \vdots\\ y_N\\ \end{bmatrix}$$
从底往上算，
$x_N*f_N=y_N$，
$x_{N-1}*f_{N-1}-c_{N-1}*x_N$
…..
$x_2*f_2-c_2*x_3$
$x_1*f_1-c_1*x_2$

计算$x_i$
则有：
$x_i = (y_i+c_i*x_{i+1})/f_i,i=N-1,N-2,...,1$

matlab代码：

function x = thomas(N, a, b, c, d)
e = zeros(N,1);
f = zeros(N,1);
x = zeros(N,1);
y = zeros(N,1);
% compute f_i and e_i
f(1) = b(1);
for i = 2:1:N
e(i) = - a(i)/f(i-1);
f(i) = b(i) + e(i)*c(i-1);
end
% compute y_i
y(1) = d(1);
for i = 2:1:N
y(i) = d(i) - e(i)*y(i-1);
end
% compute x_i
x(N) = y(N)/f(N);
for i = N-1:-1:1
x(i) = (y(i) + c(i)*x(i+1))/f(i);
end