列主元高斯消去法的实现

《数值分析课程设计》报告专业：学号：学生姓名：指导教师：一、题目列主元guess消去法求方程的解二、理论列主元高斯消去法是在高斯消去法的基础上而得到的一种比较快速合理的解线性方程组的方法。它的基本思想是每次在所在列对角线及以下元素中选择绝对值最大的元素作为主元进行消元计算。使用列主元消去法相对于高斯消去法更能减少舍入误差的影响。三、方法、算法与程序设计求解Ax=b第一步：写出增广矩阵[A|b];第二步：判断增广矩阵的秩r[A|b]与A的秩r[A]的关系：若r[A|b]=r[A]，线性方程组有唯一解；若r[A|b]r[A]，线性方程组没有解；若r[A|b]r[A]，线性方程组有无穷多解；第三步：若|A|≠0，方程组有唯一解法一：求出A-1，x=A-1b法二：利用初等行变换将[A|b]中A化为对角阵计算矩阵行列式：直接调用Det[]函数计算计算矩阵条件数：第一步：求出矩阵的逆矩阵第二步：分别计算矩阵极其逆的无穷范数，一范数和二范数第三步：求解矩阵的条件数Cond(A)∞=||A-1||∞*||A||∞Cond(A)1=||A-1||1*||A||1Cond(A)2=||A-1||2*||A||2=(λmax(A的转置*A)/λmin(A的转置*A))1/2计算机求解第一步：消元对ｋ＝１，２，３，……ｎ，进行：步骤１：选主元（第ｋ列中第ｋ个至第ｎ个元素中绝对值较大者）步骤２：将主元所在行与第ｋ行交换步骤３：消元第二步：回代求解流程图如下：k=1,2,……,ni=k+1,k+2,……,nlik=aik/akk得到aikj=k+1,k+2,……,n+1aij-aik*akj得到aij三、算例、应用实例用列主元消去法解线性方程组Ax=b⑴3.10x1+6.03x2+1.99x3=11.27x1+4.16x2-1.23x3=1;0.983x1-4.81x2+9.34x3=1⑵3.00x1+6.03x2+1.99x3=11.27x1+4.16x2-1.23x3=1.0.990x1-4.81x2+9.34x3=1分别输出A,b,detA,解向量x,⑴中A的条件数。分析比较⑴、⑵的计算结果。输出结果为：⑴A=3.1000000000000006.0300000000000001.9900000000000001.2700000000000004.160000000000000-1.2300000000000000.983000000000000-4.8100000000000009.340000000000000b=111X=-16.2344909576254896.7632694512773175.298697074088838ans=(矩阵的行列式的值)2.9967ans=（A的条件数）314.3810⑵A=3.0000000000000006.0300000000000001.9900000000000001.2700000000000004.160000000000000-1.2300000000000000.990000000000000-4.8100000000000009.340000000000000b=111X=1.0e+002*1.195273381259593-0.471426044312964-0.368402561091259ans=-0.4070虽然两题中A的数只差两个数，但是结果完全不同五、参考文献数值计算方法与算法（第二版）科学出版社数值分析（第五版）清华大学出版社六、附录%高斯列主元消元法求解线性方程组Ax=b%A为输入矩阵系数,b为方程组右端系数%方程组的解保存在x变量中formatlong;%设置为长格式显示,显示15位小数A=[？？？]det(A);cond(A);b=[？？？]'[m,n]=size(A);%先检查系数正确性ifm~=nerror('矩阵A的行数和列数必须相同');return;endifm~=size(b)error('b的大小必须和A的行数或A的列数相同');return;end%再检查方程是否存在唯一解ifrank(A)~=rank([A,b])error('A矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解');return;endc=n+1;A(:,c)=b;%(增广)fork=1:n-1[r,m]=max(abs(A(k:n,k)));%选主元m=m+k-1;%修正操作行的值if(A(m,k)~=0)if(m~=k)A([km],:)=A([mk],:);%换行endA(k+1:n,k:c)=A(k+1:n,k:c)-(A(k+1:n,k)/A(k,k))*A(k,k:c);%消去endendx=zeros(length(b),1);%回代求解x(n)=A(n,c)/A(n,n);fork=n-1:-1:1x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n))/A(k,k);enddisp('X=');disp(x);formatshort;%设置为默认格式显示,显示5位