遇到隐式格式,我们需要求解一个线性方程组。怎么办呢?当然是Thomas algorithm
注意到第一处箭头。
接着是矩阵化:
对角占有就可以LU分解
重头戏:Thomas algorithm
1.首先是计算 f1 f 1
L
的第一行乘以
U
的第一列,得到
A
的第一行第一列的元素,即
2.接着计算出每个 ei,fi,i=2,3...N e i , f i , i = 2 , 3... N
计算 e2 e 2 :
L
的第二行乘以
U
的第一列,得到
A
第二行第一列的元素,即
得到 e2=−a2/f1 e 2 = − a 2 / f 1
计算 e3 e 3 :
L
的第三行乘以
U
的第二列(乘以第一列为零),得到
A
第三行第二列的元素,即
得到 e3=−a3/f2 e 3 = − a 3 / f 2
依次类推~~~
3.计算 fi,i=2,3..N f i , i = 2 , 3.. N
L
的第二行乘以
U
的第二列,得到
A
的第二行第二列的元素,即
得到 f2=b2+e2∗c1 f 2 = b 2 + e 2 ∗ c 1
计算 f3 f 3 :
L
的第三行乘以
U
的第三列,得到
A
的第三行第三列的元素,即
得到 f3=b3+e3∗c2 f 3 = b 3 + e 3 ∗ c 2
依次类推
计算 yi,i=1,2,3,...,N y i , i = 1 , 2 , 3 , . . . , N
由
L
的
L∗Y=d L ∗ Y = d
L的第一行乘以Y,即 1∗y1=d1 1 ∗ y 1 = d 1 , y1=d1 y 1 = d 1
L的第二行乘以Y,即 e2∗y1+1∗y2=d2 e 2 ∗ y 1 + 1 ∗ y 2 = d 2 , y2=d2−e2∗y1 y 2 = d 2 − e 2 ∗ y 1
L的第三行乘以Y,即 e3∗y2+1∗y3=d3 e 3 ∗ y 2 + 1 ∗ y 3 = d 3 , y3=d2−e2∗y2 y 3 = d 2 − e 2 ∗ y 2
依次类推,
yi=di−di∗yi−1,i=2,3,...,N y i = d i − d i ∗ y i − 1 , i = 2 , 3 , . . . , N
计算
xi,i=N,N−1,...,1
x
i
,
i
=
N
,
N
−
1
,
.
.
.
,
1
由上面计算出
yi
y
i
从底往上算,
xN∗fN=yN x N ∗ f N = y N ,
xN−1∗fN−1−cN−1∗xN x N − 1 ∗ f N − 1 − c N − 1 ∗ x N
…..
x2∗f2−c2∗x3 x 2 ∗ f 2 − c 2 ∗ x 3
x1∗f1−c1∗x2 x 1 ∗ f 1 − c 1 ∗ x 2
计算
xi
x
i
则有:
xi=(yi+ci∗xi+1)/fi,i=N−1,N−2,...,1
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