|||
三对角线性方程组 $A x = b$ 在三次样条插值和微分方程边值问题的差分解法等问题中会遇到,这里讨论求解三对角方程组的追赶法的MATLAB程序。
设
$\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}\begin{bmatrix} x_1\\ x_2\\ \vdots\\ x_4\\ x_5 \end{bmatrix} = \begin{bmatrix} d_1\\ d_2\\ \vdots\\ d_{n-1}}\\ d_n \end{bmatrix}$
记为 $\textbfsymbol{A{}x} = \textbfsymbol{d}$ ,满足条件
$\left\{
\begin{array}{c}
\mid b_1\mid>\mid c_1\mid>0, \\
\mid b_i\mid\geqslant \mid a_i \mid +\mid c_i \mid, \\
\mid b_n\mid >\mid a_n\mid >0,
\end{array}
\right.
a_i{}c_i\ne0(i=2,3,\cdots,n-1).$
对 $\textbfsymbol{A}$ 进行Doolittle 分解(LU 分解)
$\begin{array}{c}
\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}\\
=\begin{bmatrix}
1&\\
l_2&1&\\
&l_3&1&\\
& &\ddots&\ddots&\\
& & &l_n &1\\
\end{bmatrix}
\begin{bmatrix}
u_1&c_1\\
&u_2&c_2\\
& &\ddots&\ddots\\
& & &u_{n-1}&c_{n-1}\\
& & & &b_n\\
\end{bmatrix}
\end{array}$
其中 $l_i,u_i$ 待定, $\textbfsymbol{U}$ 次对角线上的元素与 $\textbfsymbol{A}$ 对应次对角线上的元素相同(直接验算可知)。
利用矩阵相等可得
$\left\{ \begin{array}{c} u_1=b_1, \\ l_i = a_i/u_{i-1}, \\ u_i = b_i-l_i{}c_{i-1}, \end{array} \right. (i=2,3,\cdots,n)$
实现了 $\textbfsymbol{A}$ 的分解.
求解 $\textbfsymbol{Ly} = \textbfsymbol{d}$ , 得
$\left\{ \begin{array}{c} y_1=d_1, \\ y_i = d_i-l_i{}y_{i-1}, \end{array} \right. (i=2,3,\cdots,n)$
再求解 $\textbfsymbol{Ux} = \textbfsymbol{d}$ , 得
$\left\{
\begin{array}{c}
x_n=y_n/u_n, \\
x_i = (y_i-c_i{}x_{i+1})/u_i,
\end{array}
\right.
(i=n-1,n-2,\cdots,2,1)$
根据上面分析结果,我们用MATLAB编程可得
function [x,y,l,u] = zhuigan(a,b,c,d)
%定义函数 zhuigan
n = length(b);
%将a设置为n维向量
if length(a)==(n-1)
for i = n-1:-1:1
a(i+1)=a(i);
end
end
% 6.27 LU 分解
u = b; l = b;
for i = 2:n
l(i) = a(i)/u(i-1); u(i) = b(i) - l(i) *c(i-1);
end
% 6.28 Ly = d
y = d;
for i = 2:n
y(i) = d(i)-l(i)*y(i-1);
end
% 6.29 Ux = y
x = y; x(n) = y(n)/u(n);
for i = n-1:-1:1
x(i) = (y(i)-c(i)*x(i+1))/u(i);
end
例 用追赶法解三对角方程组
$\begin{bmatrix} 2&1\\ &1&3&1\\ & &1&1&1\\ & & &2 &1\\ \end{bmatrix} \begin{bmatrix} x_1\\x_2\\x_3\\x_4 \end{bmatrix}= \begin{bmatrix} 1\\2\\2\\0\\ \end{bmatrix}$
>> a = [1 1 2]; b = [2 3 1 1]; c = [1 1 1]; d = [1 2 2 0]; x = zhuigan(a,b,c,d)
x =
0 1.0000 -1.0000 2.0000
参考文献:
姚仰新, 王福昌, 罗家洪, 庄楚强. 高等工程数学(第三版)[M].华南理工大学出版社, 广州, 2016-08.
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-5-20 10:58
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社