追赶法求解三对角方程组

原创 2015年12月04日 09:37:50

1. 来源和背景

对于一个(主)三对角方程组,我们常用“追赶法”来进行求解. 而三对角方程组常常出现于微分方程的数值求解,例如热传导方程的边值问题

{y′′(x)=f(x,y,y), axby(a)=η1, y(b)=η2

f(x,y,y)是一个线性函数时,对该边值问题的数值解转化为一个典型的三对角方程组求解.

“追赶法”目前比较可靠的来源是下面的文章:
Thomas, L.H., Elliptic Problems in Linear Differential Equations over a Network. Watson Science Computer Laboratory Report, 1949.
其中的一个依据是,在国外的文章和教材中,“追赶法”被称为“Thomas算法”.


2. 追赶法的基本原理

追赶法的基本原理是矩阵的LU分解,即将矩阵A分解为

A=LU

其中,L为一个对角线上元素为1的下三角矩阵,U为一个上三角矩阵. 容易验证,一个三对角矩阵作LU分解以后,得到一个下二对角矩阵与一个上二对角矩阵的乘积,即
A=a11a21a12a22a32a23a33a34an1,n2an1,n1an,n1an1,nan1,n

L=1211321n1,n21n,n11

U=u11u12u22u23u33u34un1,n1un1,nun1,n

三对角矩阵A的LU分解计算过程如下:

for i = 2 to n
    A(i,i-1) = A(i,i-1)/A(i-1,i-1);
    A(i,i) = A(i,i) - A(i-1,i) * A(i,i-1);
end

在计算过程中,将下三角矩阵L和上三角矩阵U的值保存在原矩阵A中. 计算结束以后,矩阵A中的元素为

u1121u12u2232u23u33u34n1,n2un1,n1n,n1un1,nun1,n

注: 三对角矩阵A做LU分解以后,严格上三角部分的元素没有发生变化,即上三角矩阵U中的元素

ui,i+1=ai,i+1, i=1,2,,n1


3. 追赶法求解三对角方程组

使用LU分解的求解线性方程组时,不需要存储下三角矩阵,而上三角矩阵将被用于回代求解.

3.1 “追”的过程:分解

对于n阶的三对角方程组

Ax=b

我们先用LU分解得到
Ux=y

注:Ax=LUx=b,得

Ux=L1b

y=L1b,即得到方程组Ux=y.

计算过程如下:

for i = 2 to n
    A(i,i-1) = A(i,i-1)/A(i-1,i-1);
    A(i,i) = A(i,i) - A(i-1,i) * A(i,i-1);
    b(i) = b(i) - b(i-1) * A(i,i-1);
end

循环里面的前两行与LU分解完全相同,第三行负责对常数项做相应的变换. 在计算过程中,上三角矩阵U的值保存在原矩阵A中,变换后的常数y=L1b保存在b中.

3.2 “赶”的过程:回代

接着,我们用回代法求解上三角形方程组. 从三对角矩阵得到的上三角形方程组如下:

u11u12u22u23u33u34un1,n1un1,nun1,nx1x2xn1xn=y1y2yn1yn

注意在前面的计算过程中,我们将上三角矩阵U保存在A中,常数项y保存在b中. 因此,我们得到如下的回代过程:

x(n) = b(n) / A(i,i); 
for i = n-1 to 1
    x(i) = (b(i) - A(i,i+1) * x(i+1)) / A(i,i);
end

4. 实用的程序代码

在三对角矩阵中,三对角线以外的元素均为0,为了提高存储的效率,我们只需存储三对角线上的元素即可. 因此,对于前面的矩阵A,我们只存储三个向量:

d=[A(1,1),A(2,2),...,A(n,n)];
u=[A(1,2),A(2,3),...,A(n-1,n)];
l=[A(2,1),A(3,2),...,A(n,n-1)];

这三个向量分别为矩阵A三条对角线上的元素. 假定常数向量为

b=[b(1),b(2),...,b(n)];

则实用的追赶法(亦称为“Thomas算法”)求解三对角方程组的过程如下:

% 追
for i = 2 to n
    l(i-1) = l(i-1)/d(i-1);
    d(i,i) = d(i,i) - u(i-1) * l(i-1);
    b(i) = b(i) - b(i-1) * l(i-1);
end
% 赶
x(n) = b(n) / d(i); 
for i = n-1 to 1
    x(i) = (b(i) - u(i) * x(i+1)) / d(i);
end
版权声明:本文为博主原创文章,未经博主允许不得转载。

求解三对角方程组的追赶法(Matlab程序)

clear all;clc; fprintf('输入n:(10,20,30)\n'); n=input(''); n a=zeros(1,n);b=zeros(1,n);c=zeros(1,n...
  • zhangchao3322218
  • zhangchao3322218
  • 2012年03月30日 18:43
  • 10366

追赶法

给ab范围让a/b接近一个数,用追赶法。 直接枚举,取a,b为1每次对a,b求商,如果a / b > x,则a增加1,否则b增加1,每次记录下差值最小时a,b的值。...
  • ACpac
  • ACpac
  • 2014年12月30日 10:15
  • 390

追赶法

//追赶法 #include #include #define N 4 //定义一个4*4的矩阵,改变N可以改变矩阵的大小 void TDMA(float a[N],flo...
  • autumn20080101
  • autumn20080101
  • 2012年03月22日 13:39
  • 732

计算方法之用追赶法求线性方程组

/************************************* * 用追赶法求线性方程组 * * |- -| |- -| |- -| * | 2 -1...
  • wzhg0508
  • wzhg0508
  • 2013年06月11日 12:57
  • 1199

三对角阵的LU分解和三对角方程组的求解(C语言)

/*三对角阵的LU分解和三对角方程组的求解 -------------A=LU的分解算法------- 参考教材:《数值分析》李乃成,梅立泉,科学出版社     《计算方法教程》第二版 凌永...
  • zhangchao3322218
  • zhangchao3322218
  • 2012年03月30日 18:42
  • 2941

利用追赶法来求解方程Ax=b的C++程序

利用追赶法来求解矩阵Ax=b。
  • datouniao1
  • datouniao1
  • 2015年11月28日 18:56
  • 1083

数值分析 追赶法求解三对角线性方程组 MATLAB实现

函数主体部分编程算法  参考 数值分析 第四版 颜庆津 P27 运行结果截图: %追赶法求解三对角线性方程组,Ax=b,A用一维数组a,c,d存储。 function [L,U,x]=cr...
  • jingmiaa
  • jingmiaa
  • 2015年11月05日 11:48
  • 4917

最小二乘法求回归直线方程的推导过程

在数据的统计分析中,数据之间即变量x与Y之间的相关性研究非常重要,通过在直角坐标系中做散点图的方式我们会发现很多统计数据近似一条直线,它们之间或者正相关或者负相关。虽然这些数据是离散的,不是连续的,我...
  • MarsJohn
  • MarsJohn
  • 2017年02月07日 16:38
  • 21991

最小二乘法的几何意义 – 巧妙记忆公式的方法

上次写了篇文章来阐述几何投影与傅里叶级数的联系,今天我想谈谈几何投影与最小二乘法的联系,这种联系的好处是不管多复杂的公式,又可以被瞬间记住了。本文的中心思想是:最小二乘法中的几何意义是高维空间中的一个...
  • xiaoshengforever
  • xiaoshengforever
  • 2013年11月08日 20:57
  • 14536

最小二乘法拟合圆公式推导及其实现

1.1最小二乘拟合圆介绍与推导 最小二乘法(least squares analysis)是一种数学优化技术,它通过最小化误差的平方和找到一组数据的最佳函数匹配。最小二乘法是用最简的方法求得一些绝对不...
  • Jacky_Ponder
  • Jacky_Ponder
  • 2017年04月21日 15:04
  • 1822
内容举报
返回顶部
收藏助手
不良信息举报
您举报文章:追赶法求解三对角方程组
举报原因:
原因补充:

(最多只允许输入30个字)