列主元高斯消去法和列主元三角分解法解线性方程Word文档格式.docx

1、（4）回代求解1、列主元三角分解法对方程组的增广矩阵 经过k-1步分解后，可变成如下形式：第k步分解，为了避免用绝对值很小的数作除数，引进量，于是有=如果 ，则将矩阵的第t行与第k行元素互换，将（i，j）位置的新元素仍记为或，然后再做第k步分解，这时【列主元高斯消去法程序流程图】【列主元高斯消去法Matlab主程序】function x=gauss1（A,b,c） %列主元法高斯消去法解线性方程Ax=bif （length（A）=length（b） %判断输入的方程组是否有误 disp（输入方程有误！） return;enddisp（原方程为AX=b：） %显示方程组Ab-n=length（A

2、）;for k=1:n-1 %找列主元 p,q=max（abs（A（k:n,k）; %找出第k列中的最大值，其下标为p,q q=q+k-1; %q在A（k:n,k）中的行号转换为在A中的行号 if abs（p）ktemp1=A（k,: %列主元所在行不是当前行，将当前行与列主A（k,:）=A（q,: 元所在行交换（包括b） A（q,:）=temp1; temp2=b（k,: b（k,:）=b（q,: b（q,:）=temp2; %消元 for i=k+1:n m（i,k）=A（i,k）/A（k,k）; %A（k,k）将A（i,k）消为0所乘系数 A（i,k:n）=A（i,k:n）-m（i,k）

3、*A（k,k:n）; %第i行消元处理 b（i）=b（i）-m（i,k）*b（k）; %b消元处理 end消元后所得到的上三角阵是） A %显示消元后的系数矩阵b（n）=b（n）/A（n,n）; %回代求解for i=n-1:-1:1 b（i）=（b（i）-sum（A（i,i+1:n）*b（i+1:n）/A（i,i）;clear x;AX=b的解x是x=b;【调用函数解题】【列主元三角分解法程序流程图】【列主元三角分解法Matlab主程序】自己编的程序：function x=PLU（A,b,eps） %定义函数列主元三角分解法函数if （length（A）=length（b） %判断输入的方程

4、组是否有误） %显示方程组A=A b; %将A与b合并，得到增广矩阵for r=1: if r=1 for i=1: c d=max（abs（A（:,1）; %选取最大列向量，并做行交换 if c=eps %最大值小于e，主元太小，程序结束 else d=d+1-1; p=A（1,: A（1,:）=A（d,: A（d,:）=p; A（1,i）=A（1,i）; A（1,2:n）=A（1,2: A（2:n,1）=A（2:n,1）/A（1,1）; %求u（1,i） else u（r,r）=A（r,r）-A（r,1:r-1）*A（1:r-1,r）; %按照方程求取u（r,i） if abs（u（r,r

5、）=eps %如果u（r,r）小于e，则交换行 p=A（r,: A（r,:）=A（r+1,: A（r+1,: for i=r: A（r,i）=A（r,i）-A（r,1:r-1,i）; %根据公式求解，并把结果存在矩阵A中 for i=r+1:A（i,r）=（A（i,r）-A（i,1:r-1,r）/A（r,r）; %根据公式求解，并把结果存在矩阵A中 y（1）=A（1,n+1）; for i=2: h=0; for k=1:i-1 h=h+A（i,k）*y（k）; y（i）=A（i,n+1）-h; %根据公式求解y（i） x（n）=y（n）/A（n,n）; for i=n-1: for k=i+

6、1: h=h+A（i,k）*x（k）; x（i）=（y（i）-h）/A（i,i）; %根据公式求解x（i） Ax=x; %输出方程的解可直接得到P,L,U并解出方程解的的程序（查阅资料得子函数PLU1，其作用是将矩阵A分解成L乘以U的形式。PLU2为调用PLU1解题的程序，是自己编的）（）.function l,u,p=PLU1（A） %定义子函数，其功能为列主元三角分解系数矩阵A m,n=size（A）; %判断系数矩阵是否为方阵 if m=n error（矩阵不是方阵 return if det（A）=0 %判断系数矩阵能否被三角分解矩阵不能被三角分解 u=A;p=eye（m）;l=eye

7、（m）; %将系数矩阵三角分解，分别求出P,L,U m for j=i: t（j）=u（j,i）; t（j）=t（j）-u（j,k）*u（k,i）; a=i;b=abs（t（i）; for j=i+1: if babs（t（j） b=abs（t（j）; a=j; if a=i for j=1: c=u（i,j）; u（i,j）=u（a,j）; u（a,j）=c; c=p（i,j）; p（i,j）=p（a,j）; p（a,j）=c; end c=t（a）; t（a）=t（i）; t（i）=c; u（i,i）=t（i）; u（j,i）=t（j）/t（i）; u（i,j）=u（i,j）-u（i,k）

8、*u（k,j）; l=tril（u,-1）+eye（m）; u=triu（u,0）（）.function x=PLU2（A,b） %定义列主元三角分解法的函数l,u,p=PLU1（A） %调用PLU分解系数矩阵A m=length（A）; %由于A左乘p，故b也要左乘pv=b;for q=1: b（q）=sum（p（q,1:m）*v（1:m,1）; b（1）=b（1） %求解方程Ly=b for i=2:1: b（i）=（b（i）-sum（l（i,1:i-1）*b（1:i-1）;b（m）=b（m）/u（m,m）; %求解方程Ux=yfor i=m-1: b（i）=（b（i）-sum（u（i,i+1:m）*b（i+1:m）/u（i,i）;【编程疑难】这是第一次用matlab编程，对matlab的语句还不是非常熟悉，因此在编程过程中，出现了许多错误提示。并且此次编程的两种方法对矩阵的运算也比较复杂。问题主要集中在循环控制中，循环次数多了一次或者缺少了一次，导致数据错误，一些基本的编程语句在语法上也会由于生疏而产生许多问题，但是语句的错误由于系统会提示，比较容易进行修改，数据计算过程中的一些逻辑错误，比如循环变量的控制，这些系统不会提示错误，需要我们细心去发现错误，不断修正，调试。