追赶法求解三对角方程组
要求:对于给定的三对角系数矩阵和右端项,可以求
解线性代数方程组
一、 追赶法的数学理论
设系数矩阵为三对角矩阵
?b1??a2?0 A?????0??0?c1b2a3?000?c2?b3??00?000?0000?an0??0?0?? ??cn?1??bn???an?1bn?1则方程组Ax=f称为三对角方程组。
设矩阵A非奇异,A有Crout分解A=LU,其中L为下三角矩阵,U为单位上三角矩阵,记
?b1???2?0 L?????0??0?000?00???000???n?1?n0??0?0??,??0???n???1?10??01?2?001U???????000??000??????0??00?00??
???0?n?1???1??0?2?3?00?3?可先依次求出L,U中的元素后,令Ux=y,先求解下三角方程组Ly=f得出y,再求解上三角方程组Ux=y。
事实上,求解三对角方程组的2追赶法将矩阵三角分解的计算与求解两个三角方程组的计算放在一起,使算法更为紧凑。其计算公式为:
c1f1???b??,y?11?11,?1?1??对i?2,3,?,n???i?ai,?i?bi?ai?i?1,????yi?fi??iyi?1??i??xn?yn?对i?n?1,n?2,?,1???xi?yi??ixi?1?i?ci?i(*)
二、 追赶法的算法和流程图
1.预处理
生成方程组的系数ui及其除数di,事实上,按式(*)可交替生成di与ui:
d1→u1→d2→…→un?1→dn
其计算公式为
d1?b1???ui?ci/di,i?1,2,...,n?1 ?d?b?au,i?1i?1i?i?12.追的过程
顺序生成方程组右端:
yi→y2→…→yn
据式(*)的计算公式为
y1?f1/d1?i?2,3,...,n ?y?(f?ay)/d,iii?1i?i3.赶的过程
逆序得出方程组的解xi:
xn→xn?1→…→x1
其计算公式按式为
xn?yn?i?n?1,n?2,?,1 ?x?y?ux,iii?1?i三、 追赶法的Matlab实现
function x=chase(a,b,c,f)
%求解线性方程组Ax=f,其中A是三对角阵 %a是矩阵A的下对角线元素a(1)=0 %b是矩阵A的对角线元素
%c是矩阵A的上对角线元素c(N)=0 %f是方程组的右端向量 N=length(f);
x=zeros(1,N);y=zeros(1,N); d=zeros(1,N);u= zeros(1,N); %预处理 d(1)=b(1); for i=1:N-1
u(i)=c(i)/d(i);
d(i+1)=b(i+1)-a(i+1)*u(i); end