直接法解线性方程组

直接法解线性方程组

实习题目：

仿照三对角方程组的追赶法解五对角方程组，其中系数矩阵为A,右端向量为：r。将A分解为LU。其中L为下三角，U为单位上三角。A为7*7阶的矩阵，其中对角元为4 5 6 7 8 9 10。上下次三角对角线元素为1 2 3 4 5 6 ；上下第二条对角线元素为1 2 3 4 5;右端项为：1 2 3 4 5 6 7.

要求：输出系数矩阵A，右端向量r,下三角矩阵L，单位上三角矩阵U,下三角矩阵Ly=b 的解向量y,单位上三角方程组Ux=y的解（即最终的解向量。保留七位小数。

实现方法：通过MATLAB编程实现。建立MATLAB脚本文件。

首先通仿照三对角方程组的追赶法得到五对角矩阵的实现算法。

然后又MATLAB编程实现。

实验结果（MATLAB截图）：

结果分析：

通过提供的计算数据得到最终的解向量x及中间过程产生的下三角矩阵L，单位上三角矩阵U,下三角矩阵Ly=b 的解向量y。

同时为了确保算法的正确性，我还通过MATLAB的左除运算检验得使用此算法的计算结果正确。

这里由于是用MATLAB，最终结果为分数形式，考虑到精确解一般比近似解更好，因此未化成七位小数形式。

算法实现分析：

首先计算L和U的元素。由于已知L和U的特定形式（及除了对角线和上下次对角线和上下第二条对角线外，其余为0。故通过矩阵的乘法即可得到LU中元素的计算公式。（具体算法见MATLAB程序）

算法优劣点：

1.解此题时看上去要用较多的存储单元，但实际上只需存储系数矩阵A的不为0的元素。

2.A分解为LU计算完成后，后续计算x和y的“追赶过程”运算量一般来说计算量比较小。

3.此题也可用之前的LU算法求解。但此处算法与一般的LU分解的解线性方程组的算法，相比计算量小了不少。

4.对于此处特定的对称的系数矩阵A，算法还可以进一步优化。

5.由于我在此算法中A.L U的各对角值均用一个列向量表示，一个缺点在于输出A，L,U时要重新组成矩阵形式。不过优点在于减少了存储单元。

6.另一缺点是，未能将结果封装成一个文件。

后附MATLAB代码：

c=[4,5,6,7,8,9,10];d=[1,2,3,4,5,6,0];b=[0,1,2,3,4,5,6];e=[1,2,3,4,5,0,0];a=[0,0,1,2,3,4,5];

r=[1 2 3 4 5 6 7];

w=zeros(7,1);x=zeros(7,1);y=zeros(7,1);m=zeros(7,1);n=zeros(7,1);h=zeros(7,1);

w(1)=c(1);m(1)=d(1)/c(1);n(1)=e(1)/c(1);

h(2)=b(2);w(2)=c(2)-h(2)*m(1);m(2)=(d(2)-b(2)*n(1))/w(2);n(2)=e(2)/w(2);

for k=3:5

h(k)=b(k)-a(k)*m(k-2);

w(k)=c(k)-a(k)*n(k-2)-h(k)*m(k-1);

m(k)=(d(k)-h(k)*n(k-1))/w(k);

n(k)=e(k)/w(k);

end

h(6)=b(6)-a(6)*m(4);

w(6)=c(6)-a(6)*n(4)-h(6)*m(5);

m(6)=(d(6)-h(6)*n(5))/w(6);

h(7)=b(7)-a(7)*m(5);

w(7)=c(7)-a(7)*n(5)-h(7)*m(6);

y(1)=r(1)/w(1);y(2)=(r(2)-h(2)*y(1))/w(2);

for k=3:7

y(k)=(r(k)-a(k)*y(k-2)-h(k)*y(k-1))/w(k);

end

x(7)=y(7);

x(6)=y(6)-x(7)*m(6);

for k=2:6

x(7-k)=y(7-k)-m(7-k)*x(8-k)-n(7-k)*x(9-k); end

A=[c(1),d(1),e(1), 0,0,0 0;

b(2),c(2),d(2),e(2),0,0,0;

a(3),b(3),c(3),d(3),e(3),0,0;

0,a(4),b(4),c(4),d(4),e(4),0;

0,0,a(5),b(5),c(5),d(5),e(5);

0,0,0,a(6),b(6),c(6),d(6);

0,0,0,0,a(7),b(7),c(7)]

B=r'

U=[1 m(1) n(1) 0 0 0 0;

0 1 m(2) n(2) 0 0 0;

0 0 1 m(3),n(3),0,0;

0 0 0 1 m(4),n(4),0;

0 0 0 0 1 m(5),n(5);

0 0 0 0 0 1 m(6);

0 0 0 0 0 0 1]

L=[w(1),0 0 0 0 0 0;

h(2) w(2) 0 0 0 0 0;

a(3),h(3),w(3),0 0 0 0;

0 a(4) h(4) w(4) 0 0 0;

0 0 a(5) h(5) w(5) 0 0;

0 0 0 a(6) h(6) w(6) 0;

0 0 0 0 a(7) h(7) w(7)]

y

x