列主元消去法解方程组实验报告
- 格式:docx
- 大小:44.84 KB
- 文档页数:11
<数值计算方法>实验报告1.实验名称实验2 Gauss 列主元消去法2.实验题目用Gauss 列主元消去法求解线性方程组。
0.0011 2.0002 3.0003 1.0001.0001 3.7122 4.6233 2.0002.0001 1.0722 5.6433 3.000x x x x x x x x x ++=⎧⎪-++=⎨⎪-++=⎩3.实验目的加深自己对Gauss 列主元消去法的理解和认识,并且通过做实验或做练习来加强自己Gauss 列主元消去法的掌握,学会并灵活运用Gauss 列主元消去法来求解方程组。
4.基础理论-------Gauss 列主元消去法1.Gauss 列主元消去法的基本思想是:在进行第k (k=1,2,...,n-1)步消元时,从第k 列的kk a 及以下的各元素中选取绝对值最大的元素,然后通过行变换将它交换到主元素kk a 的位置上,再进行消元。
2.Gauss 列主元消去法的优点:当kk a (k=1,2,...,n-1)的绝对值很小时,用Gauss 列主元消去法来求解方程组时,可以避免所的数值结果产生较大误差或失真。
5.实验环境实验系统:Win 7实验平台:VisualC++语言6.实验过程写出算法→编写程序→计算结果Gauss 列元消去法的算法Input:方程组未知量的个数n;增广矩阵()()1,2,...,T ij A a A A An ==,其中i=1,2,…,n; j=1,2,…,n+1Output:方程组的解x1,x2,…,xn,或失败信息。
1. for i ←1ton-1 do;2. temp ←|ii a |;3. p ←I;4. for j ←i+1 to n do5. if ||ji a >temp then6. p ←j;8. end9. end10. if temp=0 then11. |return False;12. end13. if p ≠I then14. p A ⇔i A ;//i,p 两行交换15. end//列选主元16. for j ←i+1 to n do17.*j ji i A m A -ji m ←/ji ii a a ;18. j A ←*j ji i A m A -;//消元19. end7.实验结果原方程组的解为:X1=-0.490396 , x2=-0.051035 ,x3=0.3675208.附录程序清单#include<iostream.h> #include"stdio.h"#include"math.h"void main ( ){ int n=3,i,j,k,p;doubleA[10][10]={{0.001,2.000,3.000,1.000},{-1.000,3.712,4.623,2.000},{-2.0 00,1.072,5.643,3.000}},temp,m,x[100];for(i=0;i<n;i++){ //选主元temp=fabs(A[i][i]); p=i;for(k=i+1;k<n;k++)if(fabs(A[k][i])>temp){temp=fabs(A[k][i]); p=k;}if(temp==0){ printf("\n无法求解:");return;}if(p!=i)for(j=0;j<n+1;j++){ temp=A[i][j];A[i][j]=A[p][j];A[p][j]=temp;}//消元for(k=i+1;k<n;k++){ m=A[k][i]/A[i][i];for(j=i+1;j<=n;j++)A[k][j]=A[k][j]-m*A[i][j];}}//回代for(i=n-1;i>=0;i--){x[i]=A[i][n];for(j=i+1;j<n;j++)x[i]=x[i]-A[i][j]*x[j];x[i]=x[i]/A[i][i];}printf("\nx=\n");for(i=0;i<n;i++)printf("%f \n",x[i]);}。
《数值计算方法》实验报告专业:年级:学号:姓名:成绩:1.实验名称实验2高斯列主元消去法2. :用Gauss列主消去法求解线性方程组0.001*X1+2.000*X2+3.000*X3=1.000-1.000*X1+3.217*X2+4.623*X3=2.000-2.000*X1+1.072*X2+5.643*X3=3.0003.实验目的a.熟悉运用已学的数值运算方法求解线性方程—Gauss列主消去法;b.加深对计算方法技巧的认识,正确使用计算方法来求解方程;c.培养用计算机来实现科学计算和解决问题的能力。
4.基础理论列主元消去法:a.构造增广矩阵b.找到每列绝对值的最大数;c.行变换;d.消去;e.回代5.实验环境Visual C++语言6.实验过程实现算法的流程图:7.结果分析a.实验结果与理论一致;b.由于数值设置成双精度浮点型,所以初值对计算结果影响不大;c.运用程序能更好的实现计算机与科学计算的统一和协调。
8. 附录程序清单#include<stdio.h>#include<math.h>int main(){int n=3,i,j,k,p;double a[4][4];double b[4];double x[4];double m[4][4];double temp;a[1][1]=0.001; a[1][2]=2.000; a[1][3]=3.000; b[1]=1.000;a[2][1]=-1.000; a[2][2]=3.1712; a[2][3]=2.000; b[2]=2.000;a[3][1]=-2.000; a[3][2]=1.072; a[3][3]=5.643; b[3]=3.000;for(i=1;i<=n-1;i++){temp=a[i][i];p=i;for(j=i+1;j<=n;j++)if(fabs(a[j][i])>temp){temp=a[j][i];p=j;}if(temp==0)return 0;if(p!=i) //换行{for(j=1;j<=n;j++)a[0][j]=a[i][j];for(j=1;j<=n;j++)a[i][j]=a[p][j];for(j=1;j<=n;j++)a[p][j]=a[0][j];b[0]=b[i];b[i]=b[p];b[p]=b[0];}for(j=i+1;j<=n;j++){m[j][i]=a[j][i]/a[i][i];for(k=i;k<=n;k++)a[j][k]=a[j][k]-m[j][i]*a[i][k];}}if(a[n][n]==0)return 0;x[n]=b[n]/a[n][n];for(i=n-1;i>=1;i--)//回代{temp=0;for(j=i+1;j<=n;j++)temp=temp+a[i][j]*x[j];temp=b[i]-temp;x[i]=temp/a[i][i];}for(i=1;i<=n;i++)//输出结果{printf("输出结果为:x[%d]=%lf ",i,x[i]);}printf("\n");return 0;}。
实验名称:列主元消去法解方程组1 引言我们知道,高斯消去法是一个古老的解线性方程组的方法。
而在用高斯消去法解Ax=b时,其中设A为非奇异矩阵,可能出现的情况,这时必须进行带行交换的高斯消去法。
但在实际计算中即使但其绝对值很小时,用作除数,会导致中间结果矩阵元素数量级严重增长和舍入误差的扩散,使得最后的结果不可靠。
因此,小主元可能导致计算的失败,我们应该避免采用绝对值很小的主元素。
为此,我们在高斯消去法的每一步应该在系数矩阵或消元后的低阶矩阵中选取绝对值最大的元素作为主元素,保持乘数,以便减少计算过程中舍入误差对计算解的影响。
一种方式是完全主元消去法,这种消去法是在每次选主元时,选择为主元素。
这种方法是解低阶稠密矩阵方程组的有效方法,但这种方法在选取主元时要花费一定的计算机时间。
实际计算中我们常采用部分选主元的的消去法。
列主元消去法即在每次选主元时,仅依次按列选取绝对值最大的元素作为主元素,且仅交换两行,再进行消元计算。
2 实验目的和要求运用matlab编写一个.m文件,要求用列主元消去法求解方程组(实现PA=LU):要求输出以下内容:(1)计算解x;(2) L,U;(3)整形数组IP(i)(i=1,2,…,n-1)(记录主行信息)3 算法原理与流程图(1)算法原理设有线性方程组Ax=b,其中设A为非奇异矩阵。
方程组的增广矩阵为第1步(k=1):首先在A的第一列中选取绝对值最大的元素,作为第一步的主元素:,然后交换(A,b)的第1行与第i1行元素,再进行消元计算。
设列主元素消去法已经完成第1步到第k-1步的按列选主元,交换两行,消元计算得到与原方程组等价的方程组第k步计算如下:对于k=1,2,…,n-1(1)按列选主元:即确定ik使(2)如果,则A为非奇异矩阵,停止计算。
(3)如果ik≠k,则交换[A,b]第ik行与第k行元素。
(4)消元计算消元乘数满足:(5)回代求解计算解在常数项b(n)内得到。
解线性方程组的列主元素高斯消去法和LU 分解法一、实验目的:通过数值实验,从中体会解线性方程组选主元的必要性和LU 分解法的优点,以及方程组系数矩阵和右端向量的微小变化对解向量的影响。
二、实验内容:解下列两个线性方程组(1)⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎭⎫ ⎝⎛--11134.981.4987.023.116.427.199.103.601.3321x x x(2)⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----15900001.582012151526099999.23107104321x x x x 三、实验要求:(1) 用你熟悉的算法语言编写程序用列主元高斯消去法和LU 分解求解上述两个方程组,输出Ax=b 中矩阵A 及向量b, A=LU 分解的L 及U ,detA 及解向量x.(2) 将方程组(1)中系数3.01改为3.00,0.987改为0.990,用列主元高斯消去法求解变换后的方程组,输出列主元行交换次序,解向量x 及detA ,并与(1)中结果比较。
(3) 将方程组(2)中的2.099999改为2.1,5.900001改为5.9,用列主元高斯消去法求解变换后的方程组,输出解向量x 及detA ,并与(1)中的结果比较。
(4)用MATLAB的内部函数inv求出系数矩阵的逆矩阵,再输入命令x=inv(A)*b,即可求出上述各个方程组的解,并与列主元高斯消去法和LU分解法求出的解进行比较,体会选主元的方法具有良好的数值稳定性。
用MATLAB的内部函数det求出系数行列式的值,并与(1)、(2)、(3)中输出的系数行列式的值进行比较。
四、实验过程:(1)列主元高斯消去法的主程序为function [RA,RB,n,X]=liezhuY(A,b)B=[A b]; n=length(b); RA=rank(A);RB=rank(B);zhica=RB-RA;D=det(A)if zhica>0,disp('请注意:因为RA~=RB,所以此方程组无解.')returnendif RA==RBif RA==ndisp('请注意:因为RA=RB=n,所以此方程组有唯一解.')X=zeros(n,1); C=zeros(1,n+1);for p= 1:n-1[Y,j]=max(abs(B(p:n,p))); C=B(p,:);B(p,:)= B(j+p-1,:); B(j+p-1,:)=C;for k=p+1:nm= B(k,p)/ B(p,p);B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);endendb=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp endend解方程组(1)在MATLAB工作窗口输入>>A=[3.01 6.03 1.999;1.27 4.16 -1.23;0.987 -4.819.34];b=[1;1;1];[RA,RB,n,X]=liezhuY(A,b)运行后输出结果为请注意:因为RA=RB=n,所以此方程组有唯一解. D=-0.1225RA =3 RB =3 n =3X = 397.8654-157.6242-123.1120解方程组(2)在MATLAB工作窗口输入>>A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];b=[8;5.900001;5;1];[RA,RB,n,X]=liezhu(A,b)运行后输出结果为请注意:因为RA=RB=n,所以此方程组有唯一解. D=-762.0000RA =4 RB =4 n =4X =0.0000-1.00001.00001.0000LU分解法及MATLAB主程序为function hl=zhjLU(A)[n n] =size(A); RA=rank(A);D=det(A)if RA~=ndisp('请注意:因为A的n阶行列式hl等于零,所以A不能进行LU分解.A的秩RA如下:'), RA,hl=det(A);returnendif RA==nfor p=1:nh(p)=det(A(1:p, 1:p));endhl=h(1:n);for i=1:nif h(1,i)==0disp('请注意:因为A的r阶主子式等于零,所以A不能进行LU分解.A 的秩RA和各阶顺序主子式值hl依次如下:'), hl;RAreturnendendif h(1,i)~=0disp('请注意:因为A的各阶主子式都不等于零,所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下:')for j=1:nU(1,j)=A(1,j);endfor k=2:nfor i=2:nfor j=2:nL(1,1)=1;L(i,i)=1;if i>jL(1,1)=1;L(2,1)=A(2,1)/U(1,1); L(i,1)=A(i,1)/U(1,1); L(i,k)=(A(i,k)- L(i,1:k-1)*U(1:k-1,k))/U(k,k);elseU(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j);endendendendhl;RA,U,Lendend解方程组(1)在MATLAB工作窗口输入>>A=[3.01 6.03 1.999;1.27 4.16 -1.23;0.987 -4.819.34];h1=zhjLU(A)运行输出结果为请注意:因为A的各阶主子式都不等于零,所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下:D=9.8547RA =3U =3.0100 6.0300 1.99900 4.1600 -2.07340 0 5.3016L =1.0000 0 00.4219 1.0000 00.3279 -1.6316 1.0000h1 =3.0100 4.8635 -0.1225解方程组(2)在MATLAB工作窗口输入>>A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 02];h1=zhjLU(A)运行后输出结果为请注意:因为A的各阶主子式都不等于零,所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下:D=-762.0000RA =4U =10.0000 -7.0000 0 1.00000 2.1000 6.0000 2.30000 0 -2.1429 -4.23810 -0.0000 0 12.7333L =1.0000 0 0 0-0.3000 1.0000 0 00.5000 1.1905 1.0000 -0.00000.2000 1.1429 3.2000 1.0000h1 =10.0000 -0.0000 -150.0001 -762.0001(2)在MATLAB工作窗口输入>>A=[3.01 6.03 1.999;1.27 4.16 -1.23;0.987 -4.819.34];b=[1;1;1];A(1,1)=3;A(1,3)=0.990;[RA,RB,n,X]=liezhu(A,b)请注意:因为RA=RB=n,所以此方程组有唯一解.RA =3 RB =3 n =3X = -4.02641.91931.5210hi = 3.0000 4.8219 9.8547在MATLAB工作窗口输入x=[397.8654;-157.6242;-123.1120]';x1=[-4.0264;1.9193;1.5210]';wucha=x1-x运行后输出结果为wucha =-401.8918 159.5435 124.6330(3)在MATLAB工作窗口输入>>A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];A(2,2)=2.1;b(2,1)=5.9;b=[8;5.900001;5;1];[RA,RB,n,X]=lie zhu(A,b)运行后输出结果为请注意:因为RA=RB=n,所以此方程组有唯一解.RA =4 RB =4 n =4X =0.0000-1.00001.00001.0000h1 =10.0000 -0.0000 -150.0000 -762.0000在MATLAB工作窗口输入>>x=[0;-1;1;1]';x1=[0;-1;1;1]';wucha=x1-x运行后输出结果为wucha = 0 0 0 0(4)解方程组(1)在MATLAB工作窗口输入>>A=[3.01 6.03 1.999;1.27 4.16 -1.23;0.987 -4.81 9.34];B=inv(A)运行后结果为B =-268.9293 538.3418 128.4529106.7599 -213.4281 -50.956183.3992 -166.8022 -39.7090在MATLAB工作窗口输入>>b=[1;1;1];x=inv(A)*b运行后结果为x =397.8654-157.6242-123.1120在MATLAB工作窗口输入>>A=[3.01 6.03 1.999;1.27 4.16 -1.23;0.987 -4.81 9.34];A(1,1)=3;A(1,3)=0.990;B=inv(A)运行输出结果为B = 3.3424 -6.1983 -1.1705-1.3269 2.7442 0.5020-1.0365 2.0682 0.4893在MATLAB工作窗口输入>>b=[1;1;1];x=inv(A)*b运行后输出结果为x =-4.02641.91931.5210解方程组(2)在MATLAB工作窗口输入>>A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];B=inv(A) 运行后结果为B =-0.0223 -0.0984 0.1181 0.1686-0.1601 -0.1181 0.1417 0.26900.0108 0.1063 0.0724 -0.07550.1024 0.1575 -0.1890 0.1969在MATLAB工作窗口输入>>b=[8;5.900001;5;1];x=inv(A)*b运行后输出结果为x = 0-1.00001.00001.0000在MATLAB工作窗口输入>>A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];A(2,2)=2.1;B=inv(A)运行后输出结果为B =-0.0223 -0.0984 0.1181 0.1686-0.1601 -0.1181 0.1417 0.26900.0108 0.1063 0.0724 -0.07550.1024 0.1575 -0.1890 0.1969在MATLAB工作窗口输入>>b=[8;5.900001;5;1];b(2,1)=5.9;x=inv(A)*b运行后输出结果为x =-0.0000-1.00001.00001.0000五、实验结果分析:实验的数学原理很容易理解,也容易上手。
列主元消去法matlab实验报告列主元消去法是一种常用的线性方程组求解方法,它通过选取主元元素来消去其他元素,从而简化方程组的求解过程。
本文将以Matlab为工具,对列主元消去法进行实验研究,并给出相应的实验报告。
我们需要明确列主元消去法的基本原理。
列主元消去法的核心思想是选取每一列的主元素,通过消去其他元素,从而将方程组转化为上三角形或下三角形的形式。
具体来说,通过选取第一列的主元素,将第一列下方的元素消去;然后选取第二列的主元素,将第二列下方的元素消去;依此类推,直到最后一列。
这样,我们就得到了一个上(下)三角形的方程组,可以通过回代(代入)的方法求解。
接下来,我们使用Matlab编写代码,实现列主元消去法。
首先,我们需要输入一个线性方程组的系数矩阵A和常数向量b,其中A 是一个n×n的矩阵,b是一个n×1的向量。
然后,我们通过选取主元素的方式进行消去操作,得到一个上三角形的方程组。
最后,我们通过回代(代入)的方法求解方程组的解。
具体实现的代码如下所示:```matlabfunction x = gauss_elimination(A, b)n = size(A, 1); %方程组的个数% 消元过程for k = 1:n-1[~, p] = max(abs(A(k:n, k))); %选取主元素 p = p + k - 1;% 交换第k行和第p行temp = A(k, :);A(k, :) = A(p, :);A(p, :) = temp;temp = b(k);b(k) = b(p);b(p) = temp;% 消去操作for i = k+1:nfactor = A(i, k) / A(k, k);A(i, :) = A(i, :) - factor * A(k, :);b(i) = b(i) - factor * b(k);endend% 回代(代入)过程x = zeros(n, 1);x(n) = b(n) / A(n, n);for i = n-1:-1:1x(i) = (b(i) - A(i, i+1:n) * x(i+1:n)) / A(i, i);endend```接下来,我们将使用一个具体的例子来说明列主元消去法的求解过程。
数值分析上机实验报告(二)一、问题描述:利用列主元消去法求解下列方程组2X1+5X2+3X3 - 2X4=72X1- 2X2+3X3+5X4=-1X1+3X2+2X3+3X4=0X1+2X2+ X3 - 2X4=4二、算法原理:由高斯消去法知道,在消去过程中可能出现a kk(k)=0的情况,这时候消去法将无法进行,所以最好选取系数矩阵(或消元后的低阶矩阵)中绝对值最大的元素作为主元,以使高斯消去法具有较好的数值稳定性。
三、实验步骤:1、det 1;2、对于k=1,2,···,n-1(1)按列选主元|a ik.k|=max|a ik|(2)如果a i.k=0,则计算停止(det(A)=0)(3)如果i k=k则转(4)换行:a kj a ik·j(j=k,k+1,···,n)b k b ikdet -det(4)消元计算对于i=k+1,···,n○1am ik=a ik/a kk○2对于j=k+1,···,n a ij a ij—m ik*a kj○3b i b i-m ik*b ik(5)det a kk*det3、如果则计算停止(det(A)=0)4、回代求解(1)b n b n/a nn(2)对于i=n-1···,2,1bi(bi-∑aij*bj)/aii5.det ann*det四、实验框图五、源程序# include <stdio.h># include<math.h># define n 4main(){int i,j,k,l;float A[n][n],b[n],x[n],max;//输入系数矩阵及右端项for(i=0;i<n;i++)for(j=0;j<n;j++){printf("A[%d][%d]=",i,j);scanf("%f;",&A[i][j]);}for(i=0;i<n;i++){printf("b[%d]=",i);scanf("%f;",&b[i]);}//列主元消去过程for(k=0;k<n-1;k++){max=abs(A[k][k]);l=k;for(i=k+1;i<n;i++)if(abs(A[i][k])>max){max=abs(A[i][k]);l=i;} if(l>k){for(j=k;j<n;j++){max=A[k][j];A[k][j]=A[l][j];A[l][j]=max;}max=b[k];b[k]=b[l];b[l]=max;}for(i=k+1;i<n;i++){max=A[i][k]/A[k][k];for(j=k+1;j<n;j++)A[i][j]=A[i][j]-max*A[k][j];b[i]=b[i]-max*b[k];}}//回代过程x[n-1]=b[n-1]/A[n-1][n-1];for(k=1;k<n;k++){i=n-k-1;x[i]=b[i];for(j=i+1;j<n;j++)x[i]=x[i]-A[i][j]*x[j];x[i]=x[i]/A[i][i];}//输出解for(i=0;i<n;i++)printf("x[%d]=%f;",i,x[i]);getchar();}六、运行结果。
《数值分析》实验报告实验序号:实验二 实验名称: 列主元消去法解方程组 学号: 姓名:任课教师: 专业班级:)1、 实验目的:用列主元Gauss 消元法解n 阶线性代数方程组:⎪⎪⎩⎪⎪⎨⎧=+⋯++⋯⋯⋯⋯⋯=+⋯++=+⋯++nn nn 2n21n12n 2n 2221211n 1n 212111b a a a b a a a b a a a x x x x x x x x x 其基本做法是把上述方程组通过列主元Gauss 消元转化为一个等价的三角形方程组,然后再进行回代就可以求出方程组的解。
列主元消元的基本做法是选取系数矩阵的每一列中绝对值最大的作为主元,然后采取和顺序Gauss 消元法相同的步骤进行 ,求得方程组的解。
要求显示出每一个列主元以及每一大步消元后的系数矩阵),...,2,1(n k =(k )A 和常数项),...,2,1(n k =(k )b ,最后显示出方程组的解),...,2,1(n i x i =。
2、 实验内容:(1)实验分析:1. 列主元Gauss 消元法的算法思想:1. 输入增广矩阵B ;。
2. 对k =1,2,…,n ,循环:(a ) 按列选主元||:max ik ni j a a ≤≤=保存主元所在行的指标k i 。
(b ) 若a=0,则系数矩阵奇异,计算停止;否则,顺序进行。
(c ) 若k i =k 则转向(d );否则换行ki kj i b b nj a a k j k ↔=↔,...,2,1 ,(d ) 计算乘子.,...,1,/n k i a a a m ik kk ik ik +=⇒=(e ) 消元: nk i b m b b nk j i a m a a k ik i i kj ij ij ij ,...,1;:,...,1,;:+=-=+=-=3. 回代 1,...,1, ,/:1-=⎪⎪⎭⎫ ⎝⎛-=∑+=n n i a b a b b ii n i j j ij i i 用右端项b 来存放解x 。
实验报告:Gauss消元法小组成员:李岚岚、邱粉珊、缪晓浓、杨水清学号:0917020040、0917010078、0917010073、0917010112一、实验问题编写两个程序,分别利用Gauss消元法和列主元消去法求解方程组二、分析及其计算过程Gauss顺序消元法:源程序:function [x]=gaussl(A,b)[n1,n2]=size(A);n3=size(b);if n1~=n2|n2~=n3|n1~=n3disp('A的行和列的维数不同!');return;endif det(A)==0disp('系数矩阵A奇异');return;end%消元过程L=eye(n1);for j=2:n1for i=j:n1L(i,j-1)=A(i,j-1)/A(j-1,j-1);A(i,:)=A(i,:)-L(i,j-1)*A(j-1,:);b(i)=b(i)-L(i,j-1)*b(j-1);endend%回代过程x(n1)=b(n1)/A(n1,n1);for t=n1-1:-1:1for k=n1:-1:t+1b(t)=b(t)-A(t,k)*x(k);endx(t)=b(t)/A(t,t);end程序的运行以及结果:>>A=[1 2/3 1/3;9/20 1 11/20;2/3 1/3 1];>>b=[2 2 2];>> [x]=gaussl(A,b)x =1 1 1Gauss列主元消去法:源程序:function [x]=gaussll(A,b) [n1,n2]=size(A);n3=size(b);if n1~=n2|n1~=n3|n2~=n3disp('输入的方程错误!');return;endif det(A)==0disp('系数矩阵A奇异');return;endmax=zeros(n1);for m=1:n1%找主元for i=m:n1if abs(A(i,m))>maxmax=A(i,:);A(i,:)=A(m,:);A(m,:)=max;maxb=b(i);b(i)=b(m);b(m)=maxb;endend%消元过程L=eye(n1);for j=2:n1for i=j:n1L(i,j-1)=A(i,j-1)/A(j-1,j-1);A(i,:)=A(i,:)-L(i,j-1)*A(j-1,:);b(i)=b(i)-L(i,j-1)*b(j-1);endendend%回代过程x(n1)=b(n1)/A(n1,n1);for t=n1-1:-1:1for k=n1:-1:t+1b(t)=b(t)-A(t,k)*x(k);endx(t)=b(t)/A(t,t);end程序的运行以及结果:>>A=[-0.002 2 2;1 0.78125 0;3.996 5.5625 4]; >>b=[0.4 1.3816 7.4178];>>[x]= gaussll(A,b)x =1.9273 -0.6985 0.9004。
2、用列主元高斯消去法解线性方程组b =Ax .⑴ ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--34.981.4987.023.116.427.199.103.601.3⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡321x x x =⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡111 ⑵ ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--34.981.4990.023.116.427.199.103.600.3⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡321x x x =⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡111 分别输出A b A det ,,,解向量x ,(1)中A 的条件数。
分析比较(1)(2)的计算结果。
程序1(列主元高斯消去解线性方程组):#include<stdio.h>#include<math.h>#define n 3void LZYGSXQ(double a[n][n],double b[n]){double x[3],L,max,det=1;int r,t,i,j,k;for(k=0;k<n-1;k++) //选组员{{max=fabs(a[k][k]);r=k;}for(i=k+1;i<n;i++){if(fabs(a[i][k])>max)r=i;for(t=k;t<n;t++){L=a[k][t];a[k][t]=a[r][t];a[r][t]=L;}L=b[k];b[k]=b[r];b[r]=L;det=-det;}for(i=k+1;i<n;i++) //高斯消去{L=a[i][k]/a[k][k];for(j=k;j<n;j++){a[i][j]=a[i][j]-L*a[k][j];}b[i]=b[i]-L*b[k];}det=a[k][k]*det;}det=a[k][k]*det;printf("高斯消去后的方程系数\n"); //输出高斯消去后的系数矩阵for(i=0;i<n;i++){for(j=0;j<n;j++){printf("%f ",a[i][j]);}printf("%f ",b[i]);printf("\n");}printf("\n");printf("行列式的值det=%f\n",det);x[n-1]=b[n-1]/a[n-1][n-1];for(i=n-2;i>=0;i--){ //回代求解xfor(j=i+1;j<n;j++)b[i]=b[i]-a[i][j]*x[j];x[i]=b[i]/a[i][i];}printf("方程的解向量=(");for(i=0;i<=n-1;i++)printf("%f ",x[i]);printf(")\n");}void main(){double a1[3][3]={{3.01,6.03,1.99},{1.27,4.16,-1.23},{0.987,-4.81,9.34}},b1[3]={1,1,1};double a2[3][3]={{3.00,6.03,1.99},{1.27,4.16,-1.23},{0.987,-4.81,9.34}},b2[3]={1,1,1};LZYGSXQ(a1,b1);printf("\n\n\n");LZYGSXQ(a2,b2);}计算结果:方程组(1)的解向量为),,493.617725-631.911376-51592.59962(1=x ,方程组(2)的解向量为)644123.41-336653.53163940.135(2,,=x 。
列主元素消去法一、实验要求1、编程实现用列主元素消去法求解方程组二、实验目的1、进一步了解如何求解方程组2、熟悉关于C语言的一些基本编程操作三、实验内容1、用列主元素消去法求解线性方程:#include<iostream>#include<cmath>#define N 20using namespace std;void load();float a[N][N];int m;int main(){int i,j;int c,k,n,p,r;float x[N],l[N][N],s,d;cout<<"下面请输入未知数的个数m=";cin>>m;cout<<endl;cout<<"请按顺序输入增广矩阵a:"<<endl;load();for(i=0;i<m;i++){for(j=i;j<m;j++)c=(fabs(a[j][i])>fabs(a[i][i]))?j:i; /*找列最大元素*/for(n=0;n<m+1;n++){s=a[i][n]; a[i][n]=a[c][n]; a[c][n]=s;} /*将列最大数防在对角线上*/ for(p=0;p<m+1;p++)cout<<a[i][p]<<"\t";cout<<endl;for(k=i+1;k<m;k++){l[k][i]=a[k][i]/a[i][i];for(r=i;r<m+1;r++) /*化成三角阵*/a[k][r]=a[k][r]-l[k][i]*a[i][r];}}x[m-1]=a[m-1][m]/a[m-1][m-1];for(i=m-2;i>=0;i--){d=0;for(j=i+1;j<m;j++)d=d+a[i][j]*x[j];x[i]=(a[i][m]-d)/a[i][i]; /*求解*/ }cout<<"该方程组的解为:"<<endl;for(i=0;i<m;i++)cout<<"x["<<i<<"]="<<x[i]<<"\t";//system("pause");return 0;}void load(){int i,j;for(i=0;i<m;i++)for(j=0;j<m+1;j++)cin>>a[i][j];}2、LU分解法#include<stdio.h>void solve(float l[][100],float u[][100],float b[],float x[],int n) {int i,j;float t,s1,s2;float y[100];for(i=1;i<=n;i++) /* 第一次回代过程开始*/{s1=0;for(j=1;j<i;j++){t=-l[i][j];s1=s1+t*y[j];}y[i]=(b[i]+s1)/l[i][i]; }for(i=n;i>=1;i--) /* 第二次回代过程开始*/{s2=0;for(j=n;j>i;j--){t=-u[i][j];s2=s2+t*x[j];}x[i]=(y[i]+s2)/u[i][i];}}void main(){float a[100][100],l[100][100],u[100][100],x[100],b[100];int i,j,n,r,k;float s1,s2;for(i=1;i<=99;i++)/*将所有的数组置零,同时将L矩阵的对角值设为1*/ for(j=1;j<=99;j++){l[i][j]=0,u[i][j]=0;if(j==i) l[i][j]=1;}printf ("input n:\n");/*输入方程组的个数*/scanf("%d",&n);printf ("input array A:\n");/*读取原矩阵A*/for(i=1;i<=n;i++)for(j=1;j<=n;j++)scanf("%f",&a[i][j]);printf ("input array B:\n");/*读取列矩阵B*/for(i=1;i<=n;i++)scanf("%f",&b[i]);for(r=1;r<=n;r++)/*求解矩阵L和U*/{for(i=r;i<=n;i++){s1=0;for(k=1;k<=r-1;k++)s1=s1+l[r][k]*u[k][i];u[r][i]=a[r][i]-s1;}for(i=r+1;i<=n;i++){s2=0;for(k=1;k<=r-1;k++)s2=s2+l[i][k]*u[k][r];l[i][r]=(a[i][r]-s2)/u[r][r];} }。
《数值分析》实验报告实验编号:实验二课题名称:列主元Gauss消去法一、算法介绍1、输入矩阵的阶数n,方程组的增广矩阵A;2、对k=0,1,…,n-2,循环:选取列中绝对值最大的元素,将主元所在的行的元素保存在数组temp[n+1]中。
若主元为零,则系数矩阵奇异,计算停止;否则,顺序进行。
如果绝对值最大的元素就在矩阵的对角线上,则进行普通高斯消元法的第一大步,否则将方程组系数换行之后再进行普通高斯消元法的第一大步;3、然后利用回代法求解线性方程组。
二、程序代码#include<iostream>#include<cmath>#include<iomanip>using namespace std;int main(){int n=0,k=0,i=0,j=0,h=0,g=0,flag=0,i1,j1;double max=0,m=0;cout<<"***利用列主元Gauss消元法求解线性方程组***"<<endl;cout<<"请输入矩阵的阶数:"<<endl;cin>>n;double a[n][n+1];double t[n+1];double x[n];memset(a,0,sizeof(a));memset(x,0,sizeof(x));cout<<"请输入方程组的增广矩阵:"<<endl;for(i=0;i<n;i++){for(j=0;j<n+1;j++){cin>>a[i][j];}}for(k=0;k<n-1;k++){max=0;j1=0;for(i=k;i<n;i++){if(fabs(a[i][k])>max){max=fabs(a[i][k]);i1=i;j1=k;}}if(max==0){cout<<"该系数矩阵为奇异矩阵,计算停止"<<endl;flag=1;break;}else{cout<<"第"<<j1+1<<"列中绝对值最大的元素是"<<a[i1][j1]<<",在线性方程组的第"<<i1+1<<"行"<<endl;if(i1!=k){for(j=0;j<=n;j++){t[j]=a[i1][j];a[i1][j]=a[k][j];a[k][j]=t[j];}}for(i=k+1;i<=n-1;i++){m=a[i][k]/a[k][k];for(j=k;j<=n;j++)a[i][j]=a[i][j]-m*a[k][j];for(g=0;g<n;g++){for(h=0;h<n+1;h++)cout<<setiosflags(ios::fixed)<<setprecision(2)<<a[g][h]<<" ";cout<<endl;}cout<<endl;}}}if(flag==0){x[n-1]=a[n-1][n]/a[n-1][n-1] ;double sum=0;for(k=n-2;k>=0;k--){sum=0;for(i=n-1;i>=k;i--)sum+=a[k][i]*x[i];x[k]=(a[k][n]-sum)/a[k][k];}cout<<"该线性方程组的解为:"<<endl;for(i=0;i<n;i++)cout<<"x"<<i+1<<"="<<setiosflags(ios::fixed)<<setprecision(2)<<x[i]<<endl;}system("pause");return 0;}三、运算结果截屏四、算法分析列主元Gauss消元法避免了普通高斯消元法中出现的问题:遇到某个主元为零或者当主元绝对值很小时,计算将会停止或求出的结果将与其实际结果相差很远。
实验一 列主元消去法【实验内容】1.掌握列主元消去法的基本思路和迭代步骤2.并能够利用列主元的高斯消去法解任意阶数的线性方程组;3、从课后题中选一题进行验证,得出正确结果,交回实验报告与计算结果。
【实验方法与步骤】1.列主元消去法基本思路设有线性方程组Ax b =,设A 是可逆矩阵。
列主元消去法的基本思想就是通过列主元的选取将初等行变换作用于方程组的增广矩阵[]|B A b =,将其中的A 变换成一个上三角矩阵,然后求解这个三角形方程组。
2.列主元高斯消去法算法描述将方程组用增广矩阵[]()(1)|ij n n B A b a ⨯+==表示。
步骤1:消元过程,对1,2,,1k n =-L(1) 选主元,找{},1,,k i k k n ∈+L 使得,max k i k ik k i na a ≤≤= (2) 如果,0k i k a =,则矩阵A 奇异,程序结束;否则执行(3);(3) 如果k i k ≠,则交换第k 行与第k i 行对应元素位置,k kj i j a a ↔,,,1j k n =+L ;(4) 消元,对,,i k n =L ,计算/,ik ik kk l a a =对1,,1j k n =++L ,计算.ij ij ik kj a a l a =-步骤 2:回代过程:(1) 若0,nn a =则矩阵奇异,程序结束;否则执行(2);(2) ,1/;n n n nn x a a +=对1,,2,1i n =-L ,计算,11/n i i n ij j ii j i x a a x a +=+⎛⎫=- ⎪⎝⎭∑[实验程序]#include<math.h>#include<stdio.h>#include<iostream>#include<conio.h>#define NUMBER 20#define Esc 0x1b#define Enter 0x0dusing namespace std;float A[NUMBER][NUMBER+1] ,ark;int flag,n;void exchange(int r,int k);float max(int k);void message();void main(){float x[NUMBER];int r,k,i,j;char celect;void clrscr();printf("\n\nUse Gauss.");printf("\n\n1.Jie please press Enter.");printf("\n\n2.Exit press Esc.");celect=getch();if(celect==Esc)exit(0);printf("\n\n input n=");scanf("%d",&n);printf(" \n\nInput matrix A and B:");for(i=1;i<=n;i++){printf("\n\nInput a%d1--a%d%d and b%d:",i,i,n,i);for(j=1;j<=n+1;j++) scanf("%f",&A[i][j]); }for(k=1;k<=n-1;k++){ark=max(k);if(ark==0){printf("\n\nIt’s wrong!");message();}else if(flag!=k)exchange(flag,k);for(i=k+1;i<=n;i++)for(j=k+1;j<=n+1;j++)A[i][j]=A[i][j]-A[k][j]*A[i][k]/A[k][k]; }x[n]=A[n][n+1]/A[n][n];for( k=n-1;k>=1;k--){float me=0;for(j=k+1;j<=n;j++){me=me+A[k][j]*x[j];}x[k]=(A[k][n+1]-me)/A[k][k];}for(i=1;i<=n;i++){printf(" \n\nx%d=%f",i,x[i]);}message();}void exchange(int r,int k){int i;for(i=1;i<=n+1;i++)A[0][i]=A[r][i];for(i=1;i<=n+1;i++)A[r][i]=A[k][i];for(i=1;i<=n+1;i++)A[k][i]=A[0][i];}float max(int k){int i;float temp=0;for(i=k;i<=n;i++)if(fabs(A[i][k])>temp){temp=fabs(A[i][k]);flag=i;}return temp;}void message() {printf("\n\n Go on Enter ,Exit press Esc!"); switch(getch()){case Enter: main();case Esc: exit(0);default:{printf("\n\nInput error!");message();} }}【实验结果】。
数值分析实验报告之高斯列主消元法一、实验目的:清楚高斯列主元消去法与高斯主元素消去法的区别,以及它提出的必要性;掌握高斯列主元消去法的原理及推导过程,会用其解简单的线性方程组。
二、实验内容:用高斯列主元消去法解方程组⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--643.5072.1000.2623.4712.3000.1000.3000.2001.0⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡321x x x =⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡000.3000.2000.1 三、实验原理:在采用高斯消去法解方程组时,小主元可能产生麻烦,即用其做除数,会导致其他元素数量级的严重增长和舍入误差的扩散,最终使得计算的解不可靠。
故应避免采用绝对值小的主元素。
在消元之前,选择一个绝对值最大的元素作为主元,用其做除数来进行消元,这样就具有较好的数值稳定性。
这就是选主元消去法。
下面详细说明列主元素消去法。
第一步:在Ax=b 即)1()1(b x A =的系数矩阵)1(A 的第一列元素中选择一个绝对值最大的元素,不妨设为)1(1l a 。
对调)1(1j a 和)1(lj a 及)1(1b 和)1(l b (j=1,2,……,n ,1≤l ≤n)。
以)1(1l a 作为新的)1(11a 进行消元(消去对调后的第2~n 个方程中的1x )。
第k 步:(1≤k ≤n-1)设第k-1步消元过程完成,得到)()(k k b x A =,检查)(k A 中第k 列的后n-k+1个元素)(k kk a ,)(1k k k a +,…,)(k nk a ,从中选出绝对值最大者,不妨设是)(k pk a ,称它为第k 列主元素。
若p=k ,则取)(k kk a 做除数直接进行消元。
若p ≠k,则将第p 个方程与第k 个方程对调,使)(k pk a 成为新的)(k kk a ,然后以其作为除数进行消元,继续这一过程,直至得到等价的三角形方程组)()(n n b x A =,下一阶段的回代过程不变。
重庆大学学生实验报告实验课程名称计算方法开课实验室DS1421学院年级专业学生姓名学号开课时间至学年第学期1.实验目的(1)高斯列主元消去法求解线性方程组的过程(2)熟悉用迭代法求解线性方程组的过程(3)设计出相应的算法,编制相应的函数子程序2.实验内容分别用高斯列主元消去法 ,Jacobi 迭代法,Gauss--Saidel 迭代法,超松弛迭代法求解线性方程组⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡------725101391444321131243301024321x x x x 3.实验过程解:(1)高斯列主元消去法编制高斯列主元消去法的M 文件程序如下:%高斯列主元消元法求解线性方程组Ax=b%A 为输入矩阵系数,b 为方程组右端系数%方程组的解保存在x 变量中format long;%设置为长格式显示,显示15位小数A=[2,10,0,-3;-3,-4,-12,13;1,2,3,-4;4,14,9,-13]b=[10,5,-2,7]'[m,n]=size(A);%先检查系数正确性if m~=nerror('矩阵A 的行数和列数必须相同');return;endif m~=size(b)error('b 的大小必须和A 的行数或A 的列数相同');return;end%再检查方程是否存在唯一解if rank(A)~=rank([A,b])error('A 矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解');return;endc=n+1;A(:,c)=b; %(增广)for k=1:n-1[r,m]=max(abs(A(k:n,k))); %选主元m=m+k-1; %修正操作行的值if(A(m,k)~=0)if(m~=k)A([k m],:)=A([m k],:); %换行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);for k=n-1:-1:1x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n))/A(k,k);enddisp('X=');disp(x);format short;%设置为默认格式显示,显示5位运行,结果如下所示:(2)Jacobi迭代法编制迭代计算的M文件程序如下:%Jacobi迭代法求解% A为方程组的增广矩阵clc;A=[2,10,0,-3,10;-3,-4,-12,13,5;1,2,3,-4,-2;4,14,9,-13,7] MAXTIME=50;%最多进行50次迭代eps=1e-5;%迭代误差[n,m]=size(A);x=zeros(n,1);%迭代初值y=zeros(n,1);k=0;%进入迭代计算disp('迭代过程X的值情况如下:')disp('X=');while 1disp(x');for i=1:1:ns=0.0;for j=1:1:nif j~=is=s+A(i,j)*x(j);endy(i)=(A(i,n+1)-s)/A(i,i);endendfor i=1:1:nmaxeps=max(0,abs(x(i)-y(i))); %检查是否满足迭代精度要求 endif maxeps<=eps%小于迭代精度退出迭代for i=1:1:nx(i)=y(i);%将结果赋给xendreturn;endfor i=1:1:n%若不满足迭代精度要求继续进行迭代x(i)=y(i);y(i)=0.0;endk=k+1;if k>MAXTIME%超过最大迭代次数退出error('超过最大迭代次数,退出');return;endend运行该程序结果如下:(3)Gauss--Saidel迭代法编制求解程序Gauss_Seidel.m如下:%Gauss_Seidel.m% A为方程组的增广矩阵clc;format long;A=[2,10,0,-3,10;-3,-4,-12,13,5;1,2,3,-4,-2;4,14,9,-13,7][n,m]=size(A);%最多进行50次迭代Maxtime=50;%控制误差Eps=10E-5;%初始迭代值x=zeros(1,n);disp('x=');%迭代次数小于最大迭代次数,进入迭代for k=1:Maxtimedisp(x);for i=1:ns=0.0;for j=1:nif i~=js=s+A(i,j)*x(j);%计算和endendx(i)=(A(i,n+1)-s)/A(i,i);%求出此时迭代的值end%因为方程的精确解为整数,所以这里将迭代结果向整数靠近的误差作为判断迭代是否停止的条件if sum((x-floor(x)).^2)<Epsbreak;end;end;X=x;disp('迭代结果:');Xformat short;运行结果如下所示:(4)超松弛迭代法编写函数M文件如下:%SOR法求解%w=1.45%方程组系数矩阵clc;A=[2,10,0,-3;-3,-4,-12,13;1,2,3,-4;4,14,9,-13]%方程组右端系数b=[10,5,-2,7]'w=1.45;%最大迭代次数Maxtime=100;%精度要求Eps=1E-5;%以15位小数显示format long;n=length(A);k=0;%初始迭代值x=ones(n,1);y=x;disp('迭代过程:');disp('x=');while 1y=x;disp(x');%计算过程for i=1:ns=b(i);for j=1:nif j~=is=s-A(i,j)*x(j);endendif abs(A(i,i))<1E-10 | k>=Maxtimeerror('已达最大迭代次数或矩阵系数近似为0,无法进行迭代'); return;ends=s/A(i,i);x(i)=(1-w)*x(i)+w*s;endif norm(y-x,inf)<Eps%达到精度要求退出计算break;endk=k+1;enddisp('最后迭代结果:');%最后的结果X=x'%设为默认显示格式format short; 结果如下:4.实验环境及实验文件存档名实验环境:Matlab7.0文件存档名:Gauss.m,Jacobi.m,Gauss_Seidel.m,SOR.m5.实验结果及分析=1.0000=2.0000=3.0000=4.0000经过验证,高斯列主元消结果正确。
《数值分析》实验报告实验序号:实验二题目名称: 列主元Gauss消元法解n阶线性代数方程组学号: 姓名:任课教师: 马季骕专业班级:计算机科学与技术(非师范)1、实验目的:用列主元Gauss消元法解n阶线性代数方程组编写一个程序实现用列主元消元法实现解方程组的问题。
2、算法分析:其基本做法是把上述方程组通过列主元Gauss消元转化为一个等价的三角形方程组,然后再进行回代就可以求出方程组的解。
列主元消元的基本做法是选取系数矩阵的每一列中绝对值最大的作为主元,然后采取和顺序Gauss 消元法相同的步骤进行,求得方程组的解。
1. 列主元Gauss消元法的算法思想:1.输入系数矩阵A,右端项b,阶n。
2.对k=1,2,…,n,循环:(a)按列选主元保存主元所在行的指标。
(b)若a=0,则系数矩阵奇异,计算停止;否则,顺序进行。
(c)若=k则转向(d);否则换行(d)计算乘子(e)消元:3. 回代:用右端项b来存放解。
3、实验分析:建立两个数组a和b,通过循环语句将n阶增广矩阵输入进去,通过对列的循环对每一列进行消去未知数,通过n小步n大步把矩阵化简成上三角形矩阵,最后通过迭代法解得方程组得解。
3、函数分析:具体程序设计:for(i=1;i<=n;i++) //消元的第一重循环{p=0;q=0;for(m=i;m<n+1;m++){if(p<abs(a[m][i])) //比较选取列主元{p=abs(a[m][i]); //确定列主元q=m; //记录列主元所在的行序列号}}for(m=1;m<n+1;m++) //交换系数{a[0][0]=a[q][m]; //用a[0][0]做暂存单元a[q][m]=a[i][m];a[i][m]=a[0][0];}b[0]=b[q]; //交换常数项,用b[0]做暂存单元b[q]=b[i];b[i]=b[0];for(t=i;t<=n-1;t++) //具体的Gauss消元算法{w=a[t+1][i];for(j=i;j<=n;j++){a[t+1][j]=a[t+1][j]-a[i][j]*(w/a[i][i]);}b[t+1]=b[t+1]-b[i]*(w/a[i][i]);}k++;}for(i=n;i>=1;i--) //回代过程{for(j=1;j<=n;j++)v=a[i][j]*x[j]+v;x[i]=(b[i]-v)/a[i][i];v=0;}4、实验数据截频:5、程序说明:本程序在Dev C++ 环境中编译运行并且通过测试,也可以在VC++ 6.0环境中编译运行。
实验名称:列主元消去法解方程组
1引言
我们知道,高斯消去法是一个古老的解线性方程组的方法。
而在用高斯消去法解Ax=b 时,其中设A为非奇异矩阵,可能出现a kk)0的情况,这时必须进行带行交换的高斯消去法。
但在实际计算中即使a kk)0但其绝对值很小时,用a;:〉作除数,会导致中间结果矩阵A(k)元素数量级严重增长和舍入误差的扩散,使得最后的结果不可靠。
因此,小主元可能导致计算的失败,我们应该避免采用绝对值很小的主元素。
为此,我们在高斯消去法的每一步应该在系数矩阵或消元后的低阶矩阵中选取绝对值最大的元素作为主元素,保持乘数1,以便减少计算过程中舍入误差对计算解的影响。
一种方式是完全主元消去法,这种消去法是在每次选主元时,选择a(二k max a j k)0
kJ k k in」
k j n
为主元素。
这种方法是解低阶稠密矩阵方程组的有效方法,但这种方法在选取主元时要花费一定的计算机时间。
实际计算中我们常采用部分选主元的的消去法。
列主元消去法即在每次
选主元时,仅依次按列选取绝对值最大的元素作为主元素,且仅交换两行,再进行消元计算
2实验目的和要求
运用matlab编写一个.m文件,要求用列主元消去法求解方程组(实现PA=L)
111111
2 11111
3 2 1111
4 3 2 1 1 1
5 4 3 2 1 1
6 5 4 3 2 1
7 6 5 4 3 2 要求输出以下内容:
(1) 计算解x;
x.
117 1x28 1x310 1x413 1x517
1
x6
22 128 x 7
(2) L,U ;
(3) 整形数组IP (i ) (i=1,2,…,n-1 )(记录主行信息)
3算法原理与流程图
(1) 算法原理
设有线性方程组Ax=b,其中设A为非奇异矩阵。
方程组的增广矩阵为
an a12L L a1n M bi
a21a22L L a2n M b2
[代口
M
M M M
a i1
M M M M
a n1a n2L L a nn M
b n
第1步(k=1):首先在A的第一列中选取绝对值最大的元素a^,作为第一步的主元素3in max
耳!0,然后交换(A,b)的第1行与第i1行元素,再进行消元计算。
1 1 i n
设列主元素消去法已经完成第1步到第k-1步的按列选主元,交换两行,消元计算得到与原方程组等价的方程组A(k)x b(k)
J1) c ⑴
an ai2
L a(k)a1(n)M
⑵
L (2)⑵
M
a22a2k a2n b2()
O M M M [A,b] [A(k),b(k)]a kk)L a kk)M.(k)bk
a k k)u a k k)1,n M(k)
d 1
M M M M
a nk)L a(k)叫
n
M
第k步计算如下:
对于k=1,2,…,n-1
(1)按列选主元:即确定ik使
时max无0
(2) 如果a ik,k 0,则A为非奇异矩阵,停止计算
(3) 如果ik 丰k ,则交换[A , b]第ik 行与第k 行元素。
(4) 消元计算
(5) 回代求解
(3)
4程序代码及注释
for k=1:n-1
a ik
E ik
电,(i
a kk
1L ,n)
a j
b E k a
kj ,
mb.
(i ,j (i
k 1,L ,n)
k 1,L ,n)
消元乘数m ik 满足: E ik
a ik a kk
1,(i k 1,L ,n)
0 / a nn
(当 a
nn
b
(b i
j b n
n a ij b j ) / a ii ,
I
i 1
0) n 1,K ,2,1
计算解(X 1,L ,X n )在常数项 b(n)内得至叽
(2)
流程图见图1
IP(k)=k;
if i<j U(i,j)=A(i,j);
elseif i==j
L(i,j)=1; U(i,j)=A(i,j);
else
L(i,j)=A(i,j);
end
end
end
end
输入n, A , b,£
det=1.0
k=1,2,…,n-1
输出det( A) 0
按列选主元a k ,k max k i n a ik0
c0= a ik ,k
血
停机
(I
(j
c0|<
Y
ik=k
换行
a kj a ik ,j (j
k, L , n)
b k b ik,det det
消元计算
k+1,L
a ij
k+1,L
det
回代求解
b n
det
(b
a nn
a nn(当a nn0)
nn nn
n
a^ b j)/ a i ,I 1,2, L n
j i 1
* det
输岀计算解及行列式值b
(I)( I=1,2,…,n)及det
,n)
,n)
停机
,a ik
* a kj
,b i
a kk det
a ik / a kk
6
附图1
4.0000
5.0000
6.0000 0
5算例分析
1、测试示例 >> A=[ 1 2 3 4 5 6];
>> b=[3 7 11];
>> [x L U IP P]=gauss(A,b) Error! Please in put aga in! >> A=[ 1 2 3 0 0 0 4 5 6]; >> b=[3 7 11];
>> [x L U IP P]=gauss(A,b)
The equati ons have no unique soluti on! x = NaN -Inf Inf
1.0000
0.2500 1.0000
0 1.0000
0 0.7500 1.5000
IP =
3 3
P =
1 0 0
0 0 1
0 1 0
2、计算过程
(1)首先输入系数矩阵A和矩阵b
>> A=[
1 1 1 1 1 1 1
2 1 1 1 1 1 1
3 2 1 1 1 1 1
4 3 2 1 1 1 1
5 4 3 2 1 1 1
6 5 4 3 2 1 1
7 6 5 4 3 2 1];
>> b=[7 8 10 13 17 22 28];
(2)输出结果
>> [x L U IP P]=gauss(A,b)
x =
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
1.000000 0000
0.2857 1.000000 000
0.42860.8000 1.000000 00
0.57140.60000.7500 1.0000000
0.71430.40000.50000.6667 1.000000
0.1429-0.2000-0.2500-0.3333-0.5000 1.00000
0.85710.20000.25000.33330.5000 --1.0000 1.0000
7.0000 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000
0 -0.7143 -0.4286 -0.1429 0.1429 0.4286 0.7143
0 -0.8000 -0.6000 -0.4000 -0.2000 0.0000
0 -0.7500 -0.5000 -0.2500 0.0000
0 -0.6667 -0.3333 -0.0000
0 0.5000 1.0000
0 1.0000
IP =
6讨论与结论
1、时间复杂度:
>> tic;[x L U IP P]=gauss(A,b);toc
Elapsed time is 0.000856 sec on ds.
2、程序优化
初次编程时,没有考虑到给一个变量赋初值的情况。
虽然在MATLAB^变量不赋初值是完全允许的,但是由于一个变量中含有多个元素时,每次改变该数组的长度,便会增加计算机时间。
另外,给程序加上一定的判断条件及报错信息,一定程度上有程序优化的作用。
因此,本程序中的以下程序段都起到了程序优化的作用。
参考文献
[1] 易大义,沈云宝,李有法•计算方法(第2版),浙江大学出版社.p.29-53.
[2] 张琨高思超毕靖编著MATLAB2010从入门到精通电子工业出版社
0 0 0 1 0 0 0
0 0 0 0 0 1 0。