3)三对角形线性方程组
?4?100000000??x1??7???14?10000000??x??5????2????0?14?1000000??x3???13???????x00?14?1000002???4????000?14?10000??x5??6? ???????
x0000?14?1000?12???6????00000?14?100??x??14????7???x000000?14?10?4???8????0000000?14?1??x??5????9???????00000000?14?????5???x10?T?,3,0?,1,2,?3,0,1,1) x*?(2,1
二、数学原理
设系数矩阵为三对角矩阵
?b1??a2?0A????0??0?c1b2a3000c2b3000000000??0?0?? ?cn?1??bn??an?1bn?10an则方程组Ax=f称为三对角方程组。
设矩阵A非奇异,A有Crout分解A=LU,其中L为下三角矩阵,U为单位上三角矩阵,记
?b1???2?0 L????0??0?0000000??0?0??,?0???n???1?10??01?2?001U????000??000?00000??0?0?? ??n?1??1???2?300?300?n?1?n可先依次求出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,?i?ci????yfi??iyi?1?i???i?xn?yn??对i?n?1,n?2,,1??xi?yi??ixi?1三、程序设计
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
%追的过程
y(1)=f(1)/d(1); for i=2:n
y(i)=(f(i)-a(i)*y(i-1))/d(i); end
%赶的过程 x(n)=y(n); for i=n-1:-1:1
x(i)=y(i)-u(i)*x(i+1); end
?i(*)
>> a=[0,-1,-1,-1,-1,-1,-1,-1,-1,-1]; >> b=[4,4,4,4,4,4,4,4,4,4];
>> c=[-1,-1,-1,-1,-1,-1,-1,-1,-1,0]; >> f=[7,5,-13,2,6,-12,14,-4,5,-5]; >> x=chase(a,b,c,f) x =
2.0000 1.0000 -3.0000 0.0000 1.0000 -2.0000 3.0000 -0.0000 1.0000 -1.0000
四、结果分析和讨论
追赶法求解的结果为x=(2,1,-3,0,1,-2,3,0,1,-1)T。求解结果与精确解一样,这表明追赶法对于求解三对角方程组具有非常高的精度,误差非常小。算法次数也较少,不选主元也可以有效的算出精确结果,是一种计算量少而数值稳定的方法。
五、完成题目的体会与收获
通过本题的求解,深刻的理解了追赶法求解三对角方程组的算法原理。学会了追赶法的matlab编程,学会了又一种求解特殊方程组的方法。追赶法在计算量方面有着巨大的优势,因此在可能的情况下应优先使用追赶法。加深了对数值计算教材知识的理解,收获非常大。