防灾数学分享 http://blog.sciencenet.cn/u/fzmath 防灾科技学院数学教研室

博文

求解三对角线性方程组 Ax = b 的追赶法的MATLAB程序

已有 19546 次阅读 2016-10-16 21:12 |个人分类:教学辅导|系统分类:教学心得| MATLAB, 追赶法

   三对角线性方程组 $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.



https://wap.sciencenet.cn/blog-292361-1009138.html

上一篇:高等工程数学 第三版 “第七章 插值方法” 姚仰新等
下一篇:Cholesky(乔累斯基)分解与求解对阵正定方程组 Ax= b 的平方根法
收藏 IP: 27.189.55.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-5-20 14:34

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部