计算方法实验三线性方程组解法列主元高斯消去法
- 格式:docx
- 大小:155.65 KB
- 文档页数:9
实验名称:列主元消去法解方程组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)内得到。
课程设计任务书前 言回顾普通解方程组的方法,一般都是先逐个削去未知变量,最终得到只有一个未知变量的方程,解之,把得到的值回代到消去变量过程中得到的方程组,逐个求出未知变量。
这种解线性方程组的基本方法就是这里要介绍的高斯消去法。
数学上,高斯消元法(或译:高斯消去法),是线性代数中的一个算法,可用来为线性方程组求解,求出矩阵的秩,以及求出可逆方阵的逆矩阵。
当用于一个矩阵时,高斯消元法会产生出一个“行梯阵式”。
高斯消元法可以用在电脑中来解决数千条等式及未知数。
高斯消元法可以用来找出一个可逆矩阵的逆矩阵。
用关联矩阵表述网络拓扑结构,并根据厂站拓扑结构和网络拓扑结构等概念简化了电力系统的拓扑结构。
根据广义乘法和广义加法的运算规则,将改进的高斯消元算法应用于电力系统拓扑结构分析中,并引入稀疏、分块处理等技术提高了上述拓扑分析的效率。
采用上述高斯消元算法对山东电网220kV 以上的变电站进行拓扑结构分析,结果表明了运用该高斯消元法进行网络拓扑分析的正确性和有效性。
用列主元素法,选取每列的绝对值最大的元素作为消去对象并作为主元素。
然后换行使之变到主元位子上,在进行消元计算。
设)()(k k b X A ,确定第k 列主元所在位置k i ,在交换k i 行和k 行后,在进行消元,并用MATLAB 软件进行求解。
目录摘要......................................................................................... 错误!未定义书签。
第1章绪论 ........................................................................... 错误!未定义书签。
第2章高斯消元法的算法描述 (2)2.1高斯消元法的原理概述 (2)2.1.1高斯消元法的消元过程 (2)2.1.2高斯消元法的回带过程 (3)2.1.3高斯消元法的复杂度分析 (4)2.2列主高斯消元法原理简介 (5)2.2.1列主高斯消元法的消元过程 (6)2.2.2列主高斯消元法的回带过程 (6)2.2.3列主高斯消元法的算法描述 (6)第3章高斯消元法的物理应用 (9)3.1电网模型的描述 (9)3.2电网模型的问题分析 (9)3.3求解计算 (11)参考文献 (13)摘 要用列主元素高斯消去法法,选取每列的绝对值最大的元素作为消去对象并作为主元素。
解线性方程组的列主元素高斯消去法和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五、实验结果分析:实验的数学原理很容易理解,也容易上手。
全选主元高斯消去法求解实系数线性方程组实验3 全选主元高斯消去法求解实系数线性方程组全选主元的基本思想和基本原理在教材的第40页,全选主元高斯消去法求解线性代数方程组的步骤在教材的第41页。
一、算法描述参数说明:n整型变量。
方程组的阶数。
a双精度实型二维数组,体积为n*n。
存放方程组的系数矩阵,返回时将被破坏。
b双精度实型一维数组,长度为n。
存放方程组右端的常数向量;返回方程组的解向量。
本函数返回整型标志值。
若返回的标志值为0,则表示原方程组的系数矩阵奇异,输出信息“fail”;若返回的标志值不为0,则表示正常返回。
#include”stdlib.h”#include “math.h”#include “stdio.h”int rgauss(n,a,b)int n;double a[],b[];{ int *js,l,k,I,j,is,p,q;double d,t;js=malloc(n*sizeof(int)); /*开辟用于记忆列交换信息的动态空间*/l=1; /*置非奇异标志*/for (k=0;k<=n-2;k++){d=0.0;for (i=k;i<=n-1;i++) /*全选主元*/for (j=k;j<=n-1;j++){t=fabs(a[i*n+j]);if (t>0) {d=t;js[k]=j;is=i;} /*记忆行、列交换信息*/ }if (d+1.0==1.0)l=0; /*置奇异标志*/else{if (js[k]!=k)for (i=0;i<=n-1;i++) /*列交换*/{p=i*n+k;q=i*n+js[k];t=a[p];a[p]=a[q];a[q]=t;}if (is!=k){ for (j=k;j<=n-1;j++) /*行交换*/{ p=k*n+j;q=is*n+j;t=a[p];a[p]=a[q];a[q]=t;}t=b[k];b[k]=b[is];b[is]=t;}}if (l==0) /*奇异返回*/, free(js); printf(…fail\n);return(0);}d=a[k*n+k];for (j=k+1;j<=n-1;j++) /*归一化*/{ p=k*n+j; a[p]=a[p]/d;}b[k]=b[k]/d;for (i=k+1;i<=n-1;i++) /*消元*/{ for (j=k+1;j<=n-1;j++){ p=i*n+j;a[p]=a[p]-a[i*n+k]*a[k*n+j];}b[i]=b[i]-a[i*n+k]*b[k];}}d=a[(n-1)*n+n-1];if (fabs(d)+1.0==1.0) /*奇异返回*/, free(js); prin tf(“fail\n”);return(0);}b[n-1]=b[n-1]/d; /*计算解向量的最后一个分量*/ for (i=n-2;i>=0; i--) /*回代*/{ t=0.0;for (j=i+1; j<=n-1; j++)t=t+a[i*n+j]*b[j];b[i]=b[i]-t;}js[n-1]=n-1;for (k=n-1;k>=0; k--) /*恢复解向量*/if (js[k]!=k){ t=b[k]; b[k]=b[js][k]}; b[js[k]]=t;} /*解向量行交换*/ free(js); /*释放动态空间*/return(l); /*返回正常标志*/}二、用高斯消去法求解下列方程组:0.2388x0+0.2471x1+0.2568x2+1.2671x3=1.84710.1968x0+0.2071x1+1.2168x2+0.2271x3=1.74710.1581x0+1.1675x1+0.1768x2+0.1871x3=1.64711.1161x0+0.1254x1+0.1397x2+0.1490x3=1.5471主函数程序如下:/*rgaus0.c*/#include “stdio.h”#include “rgauss.c”main(){ int I;static double a[4][4]={{0.2388,0.2471,0.2568,1.2671},{0.1968,0.2071,1.2168,0.2271},{0.1581,1.1675,0.1768,0.1871},{1.1161,0.1254,0.1397,0.1490}};static double b[4]={1.8471,1.7471,1.6471,1.5471}; if (rgauss(4,a,b)!=0)for (i=0; i<=3; i++)printf(“x(%d)=%e\n”,i,b*i+);}程序运行结果为x(0)=1.04058e+00x(1)=9.87051e-01x(2)=9.35040e-01x(3)=8.81282e-01。
实验三 高斯列主元消去法一、实验目的:1、掌握高斯消去法的基本思路和迭代步骤。
2、 培养编程与上机调试能力。
二、高斯列主元消去法的基本思路与计算步骤:设有方程组Ax b =,设A 是可逆矩阵。
高斯消去法的基本思想就是僵局真的初等行变换作用于方程组的增广矩阵[]B A b = ,将其中的A 变换成一个上三角矩阵,然后求解这个三角形方程组。
列主元高斯消去法计算步骤:将方程组用增广矩阵[]()(1)ij n n B A b a ⨯+== 表示。
步骤1:消元过程,对1,2,,1k n =-(1) 选主元,找{},1,,k i k k n ∈+ 使得,max k i k ikk i n a a ≤≤= (2) 如果,0k i k a =,则矩阵A 奇异,程序结束;否则执行(3)。
(3) 如果k i k ≠,则交换第k 行与第k i 行对应元素位置,k kj i j a a ↔,,,1j k n =+ 。
(4) 消元,对,,i k n = ,计算/,ikik kk l a a =对1,,1j k n =++ ,计算 .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 =- ,计算,11/n i i n ij j ii j i x a a x a +=+⎛⎫=- ⎪⎝⎭∑三:程序流程图四:程序清单:function X=uptrbk(A,b)% A是一个n阶矩阵。
% b是一个n维向量。
% X是线性方程组AX=b的解。
[N N]=size(A);X=zeros(1,N+1);Aug=[A b];for p=1:N-1[Y,j]=max(abs(Aug(p:N,p)));%返回向量的最大值存入y,最大值的序号存入j。
C=Aug(p,:);Aug(p,:)=Aug(j+p-1,:);Aug(j+p-1,:)=C;if Aug(p,p)==0'A是奇异阵,方程无惟一解'breakendfor k=p+1:Nm=Aug(k,p)/Aug(p,p);Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1);endend% 这里用到程序函数backsub来进行回代。
西京学院数学软件实验任务书实验一实验报告一、实验名称:线性方程组高斯消去法。
二、实验目的:进一步熟悉理解Guass 消元法解法思路,提高matlab 编程能力。
三、实验要求:已知线性方程矩阵,利用软件求解线性方程组的解。
四、实验原理:消元过程:设0)0(11≠a ,令乘数)0(11)0(11/a a m i i -=,做(消去第i 个方程组的i x )操作1i m ×第1个方程+第i 个方程(i=2,3,.....n )则第i 个方程变为1)1(2)1(2...i n in i b x a x a =++ 这样消去第2,3,。
,n 个方程的变元i x 后。
原线性方程组变为:⎪⎪⎪⎩⎪⎪⎪⎨⎧=++=++=++)1()1(2)1(2)1(2)1(22)1(22)0(1)0(11)0(11... . .... ...n n nn n n n n n b x a x a b x a x a b x a x a 这样就完成了第1步消元。
回代过程:在最后的一方程中解出n x ,得:)1()1(/--=n nn n n n a b x再将n x 的值代入倒数第二个方程,解出1-n x ,依次往上反推,即可求出方程组的解:其通项为3,...1-n 2,-n k /)()1(1)1()1(=-=-+=--∑k kk n k j j k kj k k k a x a bx五、实验内容:function maintest2clcclear allA=[1 3 4;2 4 5;1 4 6];%系数矩阵 b=[1 7 6]'%常数项num=length(b)for k=1:num-1for i=k+1:numif A(k,k)~=0l=A(i,k)/A(k,k); A(i,:)=A(i,:)-A(k,:).*l; b(i)=b(i)-b(k)*l; endendendAb%回代求xx(num)=b(num)/A(num,num);for i=num-1:-1:1sum=0;for j=i+1:numsum=sum+A(i,j)*x(j);endx(i)=(b(i)-sum)/A(i,i);endxEnd六、实验结果:A =1.0000 3.0000 4.0000 0 -2.0000 -3.00000 0 0.5000b =1.00005.00007.5000x =16 -25 15。
线性方程组的几种求解方法1.高斯消元法高斯消元法是求解线性方程组的一种常用方法。
该方法的基本思想是通过对方程组进行一系列简化操作,使得方程组的解易于求得。
首先将方程组表示为增广矩阵,然后通过一系列的行变换将增广矩阵化为行简化阶梯形,最后通过回代求解出方程组的解。
2.列主元高斯消元法列主元高斯消元法是在高斯消元法的基础上进行改进的方法。
在该方法中,每次选取主元时不再仅仅选择当前列的第一个非零元素,而是从当前列中选取绝对值最大的元素作为主元。
通过选取列主元,可以避免数值稳定性问题,提高计算精度。
3.LU分解法LU分解法是一种将线性方程组的系数矩阵分解为一个下三角矩阵L 和一个上三角矩阵U的方法。
首先进行列主元高斯消元法得到行阶梯形矩阵,然后对行阶梯形矩阵进行进一步的操作,得到L和U。
最后通过回代求解出方程组的解。
4.追赶法(三角分解法)追赶法也称为三角分解法,适用于系数矩阵是对角占优的三对角矩阵的线性方程组。
追赶法是一种直接求解法,将系数矩阵分解为一个下三角矩阵L和一个上三角矩阵U,然后通过简单的代数运算即可求得方程组的解。
5.雅可比迭代法雅可比迭代法是一种迭代法,适用于对称正定矩阵的线性方程组。
该方法的基本思想是通过不断迭代求解出方程组的解。
首先将方程组表示为x=Bx+f的形式,然后通过迭代计算不断逼近x的解。
6.高斯-赛德尔迭代法高斯-赛德尔迭代法是雅可比迭代法的改进方法。
该方法在每一次迭代时,使用已经更新的解来计算新的解。
相比于雅可比迭代法,高斯-赛德尔迭代法的收敛速度更快。
7.松弛因子迭代法松弛因子迭代法是一种对高斯-赛德尔迭代法的改进方法。
该方法在每一次迭代时,通过引入松弛因子来调节新解与旧解之间的关系。
可以通过选择合适的松弛因子来加快迭代速度。
以上是一些常用的线性方程组求解方法,不同的方法适用于不同类型的线性方程组。
在实际应用中,根据问题的特点和要求选择合适的求解方法可以提高计算的效率和精度。
计算方法实验报告1课题名称用列主元高斯消去法和列主元三角分解法解线性方程目的和意义高斯消去法是一个古老的求解线性方程组的方法,但由它改进得到的选主元的高斯消去法则是目前计算机上常用的解低阶稠密矩阵方程组的有效方法;用高斯消去法解线性方程组的基本思想时用矩阵行的初等变换将系数矩阵A 约化为具有简单形式的矩阵上三角矩阵、单位矩阵等,而三角形方程组则可以直接回带求解 用高斯消去法解线性方程组b Ax =其中A ∈Rn ×n 的计算量为:乘除法运算步骤为32(1)(1)(21)(1)(1)262233n n n n n n n n n n nMD n ----+=+++=+-,加减运算步骤为(1)(21)(1)(1)(1)(25)6226n n n n n n n n n n AS -----+=++=;相比之下,传统的克莱姆法则则较为繁琐,如求解20阶线性方程组,克莱姆法则大约要19510⨯次乘法,而用高斯消去法只需要3060次乘除法;在高斯消去法运算的过程中,如果出现absAi,i 等于零或过小的情况,则会导致矩阵元素数量级严重增长和舍入误差的扩散,使得最后的计算结果不可靠,所以目前计算机上常用的解低阶稠密矩阵方程的快速有效的方法时列主元高斯消去法,从而使计算结果更加精确; 2、列主元三角分解法高斯消去法的消去过程,实质上是将A 分解为两个三角矩阵的乘积A=LU,并求解Ly=b 的过程;回带过程就是求解上三角方程组Ux=y;所以在实际的运算中,矩阵L 和U 可以直接计算出,而不需要任何中间步骤,从而在计算过程中将高斯消去法的步骤进行了进一步的简略,大大提高了运算速度,这就是三角分解法采用选主元的方式与列主元高斯消去法一样,也是为了避免除数过小,从而保证了计算的精确度计算公式1、 列主元高斯消去法设有线性方程组Ax=b,其中设A 为非奇异矩阵;方程组的增广矩阵为第1步k=1:首先在A 的第一列中选取绝对值最大的元素1l a ,作为第一步的主元素:111211212222112[,]n n n l n nn n a a a a b a a a b a a a b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦a b然后交换A,b 的第1行与第l 行元素,再进行消元计算;设列主元素消去法已经完成第1步到第k -1步的按列选主元,交换两行,消元计算得到与原方程组等价的方程组 Akx=bk第k 步计算如下:对于k=1,2,…,n -11按列选主元:即确定t 使 2如果t ≠k,则交换A,b 第t 行与第k 行元素; 3消元计算消元乘数mik 满足:4回代求解2、 列主元三角分解法 对方程组的增广矩阵 经过k -1步分解后,可变成如下形式:111max 0l i i n a a ≤≤=≠(1)(1)(1)(1)(1)1112111(2)(2)(2)(2)22222()(()1)()()()()()1,1()(,)()[,][,] k k k k nk k nk n k k k k k kk kn k k k k n k k k n nn a a a a b a a a b a a b a b b a a a +++⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥→=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦A b A b ()()max 0k k tk ik k i na a ≤≤=≠,(1,,)ik ik ik kka a m i k n a ←=-=+, (,1,,), (1,,)ij ij ik kji i ik k a a m a i j k n b b m b i k n ←+=+⎧⎨←+=+⎩⎪⎪⎩⎪⎪⎨⎧--=-←←∑+=)1,,2,1(,)(1n n i a x a b x a b x ii n i j j ij i i nnn n [,]A A b =11121,11111222,122221,11,1,1,211,11,2121,112,112,1k k k k k k k j n k k j n k k k i i i k n n kk kj kn k ik ij in i nknjk k k j k n n nnk k n a a a b A a u u u u u u y l l l l l l ll l l l u u u u u y u u u u y a a b a a b l a -------------⎡→⎣⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎦第k 步分解,为了避免用绝对值很小的数kku 作除数,引进量1111 (,1,,;1,2,,) ()/ (1,2,,;1,2,,)k kj kj km mj m k ik ik im mk kkm u a l u j k k n k n l a l u u i k k n k n -=-=⎧=-=+=⎪⎪⎨⎪=-=++=⎪⎩∑∑11(,1,,)k i ik im mk m s a l u i k k n -==-=+∑,于是有kk u =ks ;如果 ,则将矩阵的第t 行与第k 行元素互换,将i,j 位置的新元素仍记为jjl 或jja ,然后再做第k 步分解,这时列主元高斯消去法程序流程图max t ik i n s s ≤≤= ()/ 1,2,,)1 (1,2,,),kk k k t iki k ik u s s s l s s i k k n l i k k n ===++≤=++即交换前的,(且列主元高斯消去法Matlab主程序function x=gauss1A,b,c %列主元法高斯消去法解线性方程Ax=bif lengthA~=lengthb %判断输入的方程组是否有误disp'输入方程有误'return;enddisp'原方程为AX=b:' %显示方程组Abdisp'------------------------'n=lengthA;for k=1:n-1 %找列主元p,q=maxabsAk:n,k; %找出第k列中的最大值,其下标为p,qq=q+k-1; %q在Ak:n,k中的行号转换为在A中的行号if absp<cdisp'列元素太小,detA≈0';break;elseif q>ktemp1=Ak,:; %列主元所在行不是当前行,将当前行与列主Ak,:=Aq,:; 元所在行交换包括bAq,:=temp1;temp2=bk,:;bk,:=bq,:;bq,:=temp2;end%消元for i=k+1:nmi,k=Ai,k/Ak,k; %Ak,k将Ai,k消为0所乘系数Ai,k:n=Ai,k:n-mi,kAk,k:n; %第i行消元处理bi=bi-mi,kbk; %b消元处理endenddisp'消元后所得到的上三角阵是'A %显示消元后的系数矩阵bn=bn/An,n; %回代求解for i=n-1:-1:1bi=bi-sumAi,i+1:nbi+1:n/Ai,i;endclear x;disp'AX=b的解x是' x=b;调用函数解题列主元三角分解法程序流程图列主元三角分解法Matlab主程序①自己编的程序:function x=PLUA,b,eps %定义函数列主元三角分解法函数if lengthA~=lengthb %判断输入的方程组是否有误disp'输入方程有误'return;enddisp'原方程为AX=b:' %显示方程组Abdisp'------------------------'n=lengthA;A=A b; %将A与b合并,得到增广矩阵for r=1:nif r==1for i=1:nc d=maxabsA:,1; %选取最大列向量,并做行交换if c<=eps %最大值小于e,主元太小,程序结束break;elseendd=d+1-1;p=A1,:;A1,:=Ad,:;Ad,:=p;A1,i=A1,i;endA1,2:n=A1,2:n;A2:n,1=A2:n,1/A1,1; %求u1,ielseur,r=Ar,r-Ar,1:r-1A1:r-1,r; %按照方程求取ur,iif absur,r<=eps %如果ur,r小于e,则交换行p=Ar,:;Ar,:=Ar+1,:;Ar+1,:=p;elseendfor i=r:nAr,i=Ar,i-Ar,1:r-1A1:r-1,i; %根据公式求解,并把结果存在矩阵A中endfor i=r+1:nAi,r=Ai,r-Ai,1:r-1A1:r-1,r/Ar,r; %根据公式求解,并把结果存在矩阵A中endendendy1=A1,n+1;for i=2:nh=0;for k=1:i-1h=h+Ai,kyk;endyi=Ai,n+1-h; %根据公式求解yiendxn=yn/An,n;for i=n-1:-1:1h=0;for k=i+1:nh=h+Ai,kxk;endxi=yi-h/Ai,i; %根据公式求解xiendAdisp'AX=b的解x是'x=x'; %输出方程的解②可直接得到P,L,U并解出方程解的的程序查阅资料得子函数PLU1,其作用是将矩阵A分解成L乘以U的形式;PLU2为调用PLU1解题的程序,是自己编的Ⅰ.function l,u,p=PLU1A %定义子函数,其功能为列主元三角分解系数矩阵A m,n=sizeA; %判断系数矩阵是否为方阵if m~=nerror'矩阵不是方阵'returnendif detA==0 %判断系数矩阵能否被三角分解error'矩阵不能被三角分解'endu=A;p=eyem;l=eyem; %将系数矩阵三角分解,分别求出P,L,Ufor i=1:mfor j=i:mtj=uj,i;for k=1:i-1tj=tj-uj,kuk,i;endenda=i;b=absti;for j=i+1:mif b<abstjb=abstj;a=j;endendif a~=ifor j=1:mc=ui,j;ui,j=ua,j;ua,j=c;endfor j=1:mc=pi,j;pi,j=pa,j;pa,j=c;endc=ta;ta=ti;ti=c;endui,i=ti;for j=i+1:muj,i=tj/ti;endfor j=i+1:mfor k=1:i-1ui,j=ui,j-ui,kuk,j;endendendl=trilu,-1+eyem;u=triuu,0Ⅱ.function x=PLU2A,b %定义列主元三角分解法的函数l,u,p=PLU1A %调用PLU分解系数矩阵A m=lengthA; %由于A左乘p,故b也要左乘p v=b;for q=1:mbq=sumpq,1:mv1:m,1;endb1=b1 %求解方程Ly=b for i=2:1:mbi=bi-sumli,1:i-1b1:i-1;endbm=bm/um,m; %求解方程Ux=y for i=m-1:-1:1bi=bi-sumui,i+1:mbi+1:m/ui,i;endclear x;disp'AX=b的解x是' x=b;调用函数解题①②编程疑难这是第一次用matlab编程,对matlab的语句还不是非常熟悉,因此在编程过程中,出现了许多错误提示;并且此次编程的两种方法对矩阵的运算也比较复杂;问题主要集中在循环控制中,循环次数多了一次或者缺少了一次,导致数据错误,一些基本的编程语句在语法上也会由于生疏而产生许多问题,但是语句的错误由于系统会提示,比较容易进行修改,数据计算过程中的一些逻辑错误,比如循环变量的控制,这些系统不会提示错误,需要我们细心去发现错误,不断修正,调试;。
列主元高斯消去法算法流程图讲解下载温馨提示:该文档是我店铺精心编制而成,希望大家下载以后,能够帮助大家解决实际的问题。
文档下载后可定制随意修改,请根据实际需要进行相应的调整和使用,谢谢!并且,本店铺为大家提供各种各样类型的实用资料,如教育随笔、日记赏析、句子摘抄、古诗大全、经典美文、话题作文、工作总结、词语解析、文案摘录、其他资料等等,如想了解不同资料格式和写法,敬请关注!Download tips: This document is carefully compiled by theeditor.I hope that after you download them,they can help yousolve practical problems. The document can be customized andmodified after downloading,please adjust and use it according toactual needs, thank you!In addition, our shop provides you with various types ofpractical materials,such as educational essays, diaryappreciation,sentence excerpts,ancient poems,classic articles,topic composition,work summary,word parsing,copy excerpts,other materials and so on,want to know different data formats andwriting methods,please pay attention!深入理解:列主元高斯消去法的算法流程图列主元高斯消去法,也称为部分主元高斯消去法,是一种在线性代数中用于解线性方程组的数值方法。
实验报告
学院:电子信息工程
实验课程:计算方法
学生姓名:
学号:
专业班级:通信工程
实验三线性方程组解法
1 目的与要求
(1)进一步理解和掌握求线性方程组数值解的有关方法和理论。
(2)完成利用列主元高斯消去法、雅可比迭代法及高斯-赛德尔迭代法求线性方程组数值解的程序设计。
本次实验只需完成列主元高斯消去法的程序设计。
(3)比较三种算法的不同特点。
2 实验内容
通过编制程序,分别用列主元高斯消去法、雅可比迭代法及高斯-赛德尔迭代法计算如下方程组的解。
设初始值为要求满足前后两次迭代结果的差向量的
1 范数小于
3 实验原理
1)列主元高斯消去法
列主元高斯消去法就是在顺序高斯消去法的基础上,每步消元之前都要进行选主元操作,即在第k 步消元前,在第k 列的元素
中选取绝对值最大的元素,设为
,然后交换第 k 行和第 p 行,继续进行消去过程,直到获得上三角方程组,然后通过回代得到方程的根。
4 程序设计
(1)流程图
列主元高斯消去程序流程图
(2)程序代码
#include<stdio.h>
#include<math.h>
void main()
{
float a[3][4],x,s;
int i,j,m,k;
printf("please input coeffient martix array:\n");
for(i=0;i<3;i++) //输入增广矩阵// {
for(j=0;j<4;j++)
{
scanf("%f",&a[i][j]);
}
}
printf("\n");
printf("Output the input matrix");
printf("\n");
for(i=0;i<3;i++) //输出输入的矩阵// {
for(j=0;j<4;j++)
{
printf("%8.4f",a[i][j]);
}
printf("\n");
}
printf("\n");
for(k=0;k<=2;k++) //在不同列中选主元// {
m=k;
for(i=k+1;i<=2;i++)
{
x=fabs(a[m][k]);
if(a[i][k]>x)
{
m=i;
}
}
if (a[m][k]==0)
{
printf("ERROR");
return;
}
else
{
for(j=k;j<=3;j++) //交换两行// {
x=a[k][j];
a[k][j]=a[m][j];
a[m][j]=x;
}
for(i=k+1;i<=2;i++) //消元过程//
{
s=a[i][k]/a[k][k];
for(j=0;j<=3;j++)
{
a[i][j]=a[i][j]-s*a[k][j];
}
}
}
}
printf("output matrix after gauss transform");
printf("\n");
for(i=0;i<3;i++) //输出变换后的矩阵//
{
for(j=0;j<4;j++)
{
printf("%8.4f",a[i][j]);
}
printf("\n");
}
printf("\n");
a[2][3]=a[2][3]/a[2][2]; //回代过程//
for(i=1;i>=0;i--)
{
x=0;
for(j=i+1;j<3;j++)
{
x=x+a[i][j]*a[j][3];
}
a[i][3]=(a[i][3]-x)/a[i][i];
}
printf("output roots of linear equations"); printf("\n");
printf("x1=%f\n",a[0][3]); //输出//
printf("x2=%f\n",a[1][3]);
printf("x3=%f\n",a[2][3]);
}
5 实验结果与分析
分析:
(1)本次试验采用了for循环及if-else判断语句,较简单的实现了线性方程组列主元高斯消去求根的编程任务。
(2)由本次试验结果来看,列主元高斯消去法避免了选用小元或零元作分母,数值稳定性较好。
对于雅克比迭代法,其计算公式简单,
每迭代一次只需计算一次矩阵和向量的乘法,且计算过程中原始矩阵A始终不变,比较容易并行计算。
然而雅克比迭代法收敛速度较慢,而且占据的存储空间较大,所以工程中一般不直接用雅克比迭代法,而用其改进方法。
对于高斯赛德尔迭代法,达到相同精度时比雅克比迭代法收敛快。
(3)此次试验较好的完成了任务,巩固了课堂学习知识。