埋深 标高 上覆岩建筑场地覆盖层厚度度之间是什么关系

线性方程组的直接解法
线性方程组的直接解法
范文一:云南大学数学与统计学实验教学中心实验报告一、实验目的练习线性方程组的直接解法。二、实验内容求解如下三对角线方程组:21??x1??2??1??2?2?1??x????3????2?????2?3?2????x3?????1??三、实验环境PC计算机、C语言四.实验方法列主元消元法五、实验过程1实验步骤(1)编程: 根据所用算法及选用语言编出源程序(2). 开机, 打开所用语言系统输入所编源程序.(3). 调试程序, 修改错误直至能正确运行.(4). 运行程序并输出计算结果.2 关键代码及其解释#define M 3
#define N 3
/*宏定义行列式有3行3列a[M][N],b[M],m[M][N],x[N]
/*分别用于存放行列式,bi、乘数及x的值; int i,j,k,l,H
/*用于存放中间变量的值3 调试过程通过动态调试f7和ctrl+f7来分步查看结果;或者通过ctrl+f7和alt+f5编译链接可得出结果,得出结果与计算一致。六、实验总结1.遇到的问题及解决过程常忘记语句结束要加;,通过编译发现错误,并改正;函数必须先定义后使用;编译链接成功后,得出结果与自己计算结果一致。2.产生的错误及原因分析Compiling D:\TC20\NONAME.C:Error D:\TC20\NONAME.C 12: S in function mainWarning D:\TC20\NONAME.C 16: Code has no effect in function main以上以上为语法错误,经过数遍的调试后就可以了!3.体会和收获。学习使用线性方程组的直接解法,即高斯消元法,同时使用编程来求解使得求解根的问题变简单,也避免了人工操作计算量大的问题以及可以避免方程组中方程个数很多很难计算的问题。高斯-若尔当消元法是高斯消元法的一种变形与改进,此算法不用换行、换列,不用回代,精度高,但缺点是循环语句比较难组织;其运算量超过高斯消元法,所以一般不用其解线性方程组,但用它求逆矩阵还是比较方便的。七:程序源代码:#include#include#define M 3#define N 3int main(){ int a[M][N],b[M],m[M][N],x[N];int i,j,k,l,H;printf("input ge xiang xi shu:\n");for(i=0;ifor(j=0;jscanf("%f\t",&a[i][j]);printf("input b[i]:\n");for(i=0;iscanf("%f\t",&b[i]);if(a[0][0]!=0)for(k==0;kfor(i=k+1;i{m[i][0]=a[i][0]/a[0][0]; for(j==0;j{a[i][j]=a[i][j]-m[i][0]*a[i][j]; b[i]=b[i]-m[i][0]*b[0];} }x[N-1]=b[N-1]/a[N-1][N-1]; printf("x[3]=:%f\n",x[N-1]); for(l=1;i>0;l++){i=N-1-l;H=a[i][i+1]*x[i];x[i]=(b[i]-H)/a[i][i];printf("x[%d]=%f\n",i,x[i]); }}运行结果:x[1]=1.00000x[2]=-1.00000x[3]=3.000000八.教师评语:原文地址:云南大学数学与统计学实验教学中心实验报告一、实验目的练习线性方程组的直接解法。二、实验内容求解如下三对角线方程组:21??x1??2??1??2?2?1??x????3????2?????2?3?2????x3?????1??三、实验环境PC计算机、C语言四.实验方法列主元消元法五、实验过程1实验步骤(1)编程: 根据所用算法及选用语言编出源程序(2). 开机, 打开所用语言系统输入所编源程序.(3). 调试程序, 修改错误直至能正确运行.(4). 运行程序并输出计算结果.2 关键代码及其解释#define M 3
#define N 3
/*宏定义行列式有3行3列a[M][N],b[M],m[M][N],x[N]
/*分别用于存放行列式,bi、乘数及x的值; int i,j,k,l,H
/*用于存放中间变量的值3 调试过程通过动态调试f7和ctrl+f7来分步查看结果;或者通过ctrl+f7和alt+f5编译链接可得出结果,得出结果与计算一致。六、实验总结1.遇到的问题及解决过程常忘记语句结束要加;,通过编译发现错误,并改正;函数必须先定义后使用;编译链接成功后,得出结果与自己计算结果一致。2.产生的错误及原因分析Compiling D:\TC20\NONAME.C:Error D:\TC20\NONAME.C 12: S in function mainWarning D:\TC20\NONAME.C 16: Code has no effect in function main以上以上为语法错误,经过数遍的调试后就可以了!3.体会和收获。学习使用线性方程组的直接解法,即高斯消元法,同时使用编程来求解使得求解根的问题变简单,也避免了人工操作计算量大的问题以及可以避免方程组中方程个数很多很难计算的问题。高斯-若尔当消元法是高斯消元法的一种变形与改进,此算法不用换行、换列,不用回代,精度高,但缺点是循环语句比较难组织;其运算量超过高斯消元法,所以一般不用其解线性方程组,但用它求逆矩阵还是比较方便的。七:程序源代码:#include#include#define M 3#define N 3int main(){ int a[M][N],b[M],m[M][N],x[N];int i,j,k,l,H;printf("input ge xiang xi shu:\n");for(i=0;ifor(j=0;jscanf("%f\t",&a[i][j]);printf("input b[i]:\n");for(i=0;iscanf("%f\t",&b[i]);if(a[0][0]!=0)for(k==0;kfor(i=k+1;i{m[i][0]=a[i][0]/a[0][0]; for(j==0;j{a[i][j]=a[i][j]-m[i][0]*a[i][j]; b[i]=b[i]-m[i][0]*b[0];} }x[N-1]=b[N-1]/a[N-1][N-1]; printf("x[3]=:%f\n",x[N-1]); for(l=1;i>0;l++){i=N-1-l;H=a[i][i+1]*x[i];x[i]=(b[i]-H)/a[i][i];printf("x[%d]=%f\n",i,x[i]); }}运行结果:x[1]=1.00000x[2]=-1.00000x[3]=3.000000八.教师评语:
范文二:第五讲 线性方程组的直接接法给定一个线性方程组Ax=b其中,?a11?aA=?21?#??an1a12?????a22,x=??,b=?2??#??#?#%#??????an2A为系数列向量,非奇异,b为右端向量,x为需求解的未知向量。解线性方程组数值解法有两类:直接法:按求精确解的方法运算求解,有Gauss消去法及修正(矩阵分解法)等。 迭代法:给一个初始近似解,按一定法则逐步求更精确的近似解的过程;有经典与现代迭代法。5.1.Gauss消去法 (Elimination Method)
2.1.1.三角形方程组的解法三角形方程组是最容易求解的,而Gauss消去法是把一般线性方程组化成两个三角形方程组来求解的。现在考虑上三角形方程组u1n??x1??b1??u11u12???????%##??#?=?#? ??????uuxbn-1,n-1n-1,n??n-1???n-1?????unn?????xn??bn?其中,uii≠0,i=1,2,???,n。于是,bxn=n,xi=unnbi-j=i+1∑uxijnjuii,i=n-1,n-2???,1同理,对于下三角形的方程组?l11??x1??b1???????llxb22?21??2??2??#??#?=?#? #%??????ullxbn-1,n-1?n-1,1n-1,2??n-1??n-1??????uln2ln,n-1lnn??n1??xn??bn?其中,lii≠0,i=1,2,???,n。于是,bx1=1,xi=l112.1.2.Gauss消去法初始增广矩阵为bi-∑lijxij=1i-1lii,i=2,???,n?a11?a(A,b)=?21?#??an1a12b1??b2? #??bn?第一步消元过程:不妨假设a11≠0(如若不然,则在该列中选个非零元,将该非零元所在行与第一行对调。) 把第1列第2到第n个元素变成0,相当于左乘行Gauss矩阵?1??-a21?aL1(l1)=?11?#?an1?-?a11?a11??0?(A,b)→??#??0?1aij=aij-a12aa22-21a12a11%an2-an1a12???1?? %??1??a1nb1??a21a21??a11a2n-a1nb2-b1a11a11??0?阅读详情:??##??#?0aaann-n1a1nbn-n1b1??a11a11?a12a122b1?1?b2?#?1?bn?ai1a1j,i,j=2,3,???,n a11第二步消元过程:不妨假设a122≠0(如若不然,则在该列主对角元以下选个非零元,将该非零元所在行与第二行对调。) 把第2列第3到第n个元素变成0,相当于左乘行Gauss矩阵。?1?1?1?a32-1?a22L2(l2)=??#?a1?2-n?1a22?????%?? 1???1???a11a12?1a22???a11a12?????a11a12?1a22??阅读详情:??????2ij1ija13a123133a1,na12n11a32a3211a-1a23a22a22a1aa-1a23a22a22143#1n3%#a1a1121n21a-1a231a22a22a13a#2an3???1a321?1b3-1b2?a22??a11421b4-1b2?a22??#?1an21?1bn-1b2?a22?b11b2a1,na12n2a3n2a4n%#2b1?1?b2?2b3?2?b4?#??2?bn?ai121其中,a=a-1a2j,i,j=3,4,???,n。a22一直进行下去,直到将A变换为一个上三角矩阵。之后运用上三角方程组的数值解法即可!2.1.3.列主元Gauss消去法在消元过程中,常出现主对角元绝对值较小或为0的情况,克服这一困难的办法是列主元消去法。列主元消去法的思想:每次消元过程先在当前变换的列元素中选绝对值最大的为主元,并根据需要交换相关的行,然后再消元。 例1.解方程组?12x1-3x2+3x3=15??18x1-3x2+x3
=15 ? -x+2x+x
=623?1【解答】 增广矩阵为?12-3315???(A,b)=?18-3115???-1216??法一 Gauss消去法??15????12--33?r2-2r1??37?7?12r3(A,b)=?18-r?2?????-1216?→??r103+12r1??222→?0?29??075444?????0因此,x3=16=3, 3-15+7x-15+73?3x2===2,22x15-3x3+3x215-3?3+3?21=12=12=1?方程组的解为x=?1??2??。??3??法二 列主元Gauss消去法-33372-2016315??-15??2?16?????18-312?12-3315?rr?18-?7??12??3(A,b)=?18-3115?→?12-3315?→?0-11?3?-1216??-r1?????1119?0618???18-31r2r3?1119
→?0?618?7?0-13???15?6?18-3?r+r?41?3112?11→0?6?6???005???119183211?15??41?6?96??11??15??5??41??6?因此,96x3==3,-x3-?3x2===2,11116615-x3+3x215-3+3?2x1===1,1818?1???因此,方恒组的解为x=?2?。?3???5.2.矩阵分解法(三角分解法)求解线线方程组Ax=b,其中A可逆。Crout分解,LDU三角分解法的思想是将A进行LU分解(包括Doolitel分解,分解,LLT分解,LDLT分解等)。设A=LU,Ax=b,先解方程组Ly=b,再解方程组Ux=y。设A=LDU,Ax=b,先解方程组Lz=b,再解方程组Dy=z,再解方程组Ux=y。若A为正定矩阵,则可以进行LLT分解和LDLT分解。A=LLT,先解Ly=b,再解LTx=y。A=LDLT,先解Lz=b,再解Dy=z,再解LTx=y。例2.求解线性方程组?3??2?1??46??x1??9??????6??x2??7?=238??x3??1??????9412??x4??13?6352【解答】?36361000?????????30102-00?3??0??12380010???→??0026-1010???→??494120001???3????0??0104-43001??????因此,?1000??3636????2?100636??2526???3102???1238??3??0?=??1010??026?? ?49412???3??0???002???4??0?3101??解方程组??00???y?1??9?y???7???1??23010??y=??3??1?????4??y????4??3101?13?????y1????9?得到?y2??=?1??。?y3??y???-2??0?4??再求解方程组636-231-130-43-100??00??10???01????3??0?0??0636??x1??9??????102??x2??1?=?????-2026x3?????002??x4??0???x1????2????x1????2??得到?x2??x?=?1?。这就求得方程组的解为?x2?=?1?。3??x???-1????x3????-1??4??0??x4??0?例3.求解线性方程组??????x1??1??x???2???2122915?=?13? ?191536????x3??41??x???4??72?【解答】??1321?
?????????040?→????1?6???1????1000??0466-3100????1?0466-3100? ?→??00164-2-310??→?0164-2-3?2???0210????00426-1-3201???????0025-1-32-1?41????1000??1000??13??1321??300????1321??00?000???66??31??1400???10???0164??=??23010????0??01?2122915??=??23?00160??91534??2?1???1??0025??2???31???00025??00???13241???0??1241???0???00??1000??1??3100????100??100?1320?3?000??1321?=?10?200??200????12??200???23??004??004??01320??=?31?34??0233??20??0???1??0??0??0041???1????000501??0005??0???4??2??315??0005???1324??00??11?21?33?22??11?4??01????1000??3100???z1??1?z1??1?解方程??2310????z????????2?=?13??2??31???z3???,得到?z2?=?10?。?41????z3???24??121???z4??72??z???50?4?4??1???1000??y1??1?解方程?400???00160???y1??5???y???????2?y?=?10?,得到?y2?=?2?3??24??y?3?。 3??00025????y???????4??50??y4???2??2????1321??03???1??13?x1解方程?22????x1??x??52???????3?-2??2?=?,得到?x2?=??1??x3??x?1??。这就求得方程组的解为3??001?4????0?x??3??2??????4??001??x4????2??2????x1??3??x???2??=?x??-2?1??。 3?x???4??2?练习题一.用高斯消去法解线性方程组?12-33??x1??1.??18-31?????x??15??2?=?-121???15? ??x3????6???140??2.??152??x1???x???2??2?=?3? ??011????x3????1????121-2?3.?253-2??x1??4?x2?????????-2-235?x??7?=?-1?? ?132-3???3?x???4??0???1103??x1??4?4.?21-11?????x???2??3-1-13??x=?1? 3??-3?-123-1????x????4??4?二.用Crout分解法解线性方程组?121??1.??223??x1??0???x??3??2?-1-30????x?=3?????2? ???223??x1?2. ??477?????x???3??2?-245????x?=3??1????-7????3636??x1??93.?2526??238????x????2??7??1?49412????x3?=?1? ?x???4??13?三.用Cholesky分解法解线性方程组?121??x1?1.??286?????x???6??2?1623????x?=3??30??? ?60????1324??x1??13?2.?313818??x??????2?=?53? ?81930????x3??41??x???4??78?四.用LDLT分解法解线性方程组?121??x1??1. ??286????x??6??2?1623????x?=3??30???60??????1324??x2. ?313818???1??13??x???2?=?53??282119?? ?4181930????x3??x???41??4??78?
范文三:求解线性方程组的直接解法5.2
Gauss消去法实现了LU分解顺序消元结束时的上三角矩阵U和所用的乘数,严格下三角矩阵。将下三角矩阵的对角元改成1,记为L,则有A=LU,?2?4????22743??1??7?2??5?????112??2????1????233??1?6??这事实是一般的,我们不难从消去的第k个元素时的矩阵k行及k列元素的历史得到这一点.因为从消元的历史有
ukj=akj-mk1u1j- mk2u2j -…- mk,k-1uk-1,j, j=k,k+1,…,n
mik=(aik-mi1u1k- mi2u2k -…-mi,k-1uk-1,k)/ukk i=k+1,k+2,…,n 于是
akj=mk1u1j+mk2u2j+…+mk,k-1uk-1,j+ukj, j=k,k+1,…,n
aik=mi1u1k+mi2u2k+…+mi,k-1uk-1,k+mikukk i=k+1,k+2,…,n 从前面两个式子我们可以直接计算L和U(见下段).将矩阵分解为单位下三角矩阵和上三角矩阵之积称为矩阵的LU分解.顺序消元实现了LU分解,同时还求出了g, Lg=b的解.②
直接LU分解上段我们得到(lij=mij)
ukj=akj-lk1u1j-lk2u2j -…- lk,k-1uk-1,j, j=k,k+1,…,n
lik=(aik-li1u1k-li2u2k -…-li,k-1uk-1,k)/ukk i=k+1,k+2,…,n2诸元素对应乘积,只不过算L的元素时还要除以同列对角元.这一规律很容易记住.可写成算法(L和U可存放于A):
for k=1:n-1
ukj=akj-lk1u1j-lk2u2j -…- lk,k-1uk-1,jendfor i=k+1:nlik=(aik-li1u1k-li2u2k -…-li,k-1uk-1,k)/ukk
end这一算法也叫Gauss消去法的紧凑格式,可一次算得L,U的元素,不需逐步计算存储.考察上面的表格会发现还可安排其它计算次序,只要在这一次序下每个元素左边的L的元素与上方的U的元素已计算在先。例如, 逐行自左而右的次序, 逐列自上而下的次序, … 易知g的计算规律同U.利用LU分解解Ax=b分三步:1.分解A=LU 2.解 Lg=b 求g 3.解 Ux=y 求x例3.
用直接LU分解法解?1??2?3?2513??x1??2??x2?5???x3??14???????18? ??20????解
用分解公式计算得?1?A??2?3?01?50??1??0??0?1???02103???4??LU. ?24??求解Ly?(14,18,20),得 y?(14,?10,?72),Ux?(14,?10,?72),得 x?(1,2,3).TTTT③
其它分解我们用顺序消元和直接分解两种方法实现了LU分解.还有更一般的三角分解,比如,下三角矩阵和单位上三角矩阵之积,又如单位下三角矩阵,对角矩阵,单位上三角矩阵之积,等等.下面给出第二种分解形式的算法LDR分解法。A=LDR,L是单位下三角矩阵,D是对角矩阵,R是单位上三角矩阵.逐列计算(逐列作LU分解,再用U的对角元素除各行),结果存入A。for j=1:n
for i=2:jaij=aij-ai1a1j-ai2a2j -…-ai,i-1ai-1,jendfor i=j+1:naij=(aij-ai1a1j-ai2a2j -…-ai,j-1aj-1,j)/ajj
for i=1:j-1
aij= aij/aii
列主元素的LU分解对照顺序消元和LU分解,列主元素法也可得列主元素的LU分解:
PA=LUP是行交换结果的排列阵,L和U同前.例4.
列主元素法解方程组并写出系数矩阵的LU分解.?2??1?2?211021?26???9???(?1/2)?(1)7???20?1021?26???6???(1)?(?1/2)1???2?100126??1????1??2???6??3???括号内是乘数,k=2时2,3行交换.因而有?1??????2??11??????22110??1??2?1??1?????1/210??2????1????2?10??1?2??直接作列主元素LU分解,因为在k步要先选主元素,所以作如下改变:
for k=1:n-1for i=k:naik=aik-li1u1k-li2u2k -…-li,k-1uk-1,k
找p:apk?maxakkak?1,k,?,ank1??endp行?k行 ik=pfor j=k+1:n
ukj=akj-lk1u1j-lk2u2j -…- lk,k-1uk-1,j endfor i=k+1:n
lik= aik/ukk end可将lik存于aik,ukj存于akj.二
实验部分 本章实验内容:实验题目:Gauss消元法,追赶法,范数。 实验内容:①编制用Gauss消元法求解线性方程组Ax=f的程序。②编制用追赶法求解线性方程组Ax=f的程序。③编制向量和矩阵的范数程序。实验目的:①了解Gauss消元法原理及实现条件,熟练掌握Gauss消元法解方程组的算法,并能计算行列式的值。②掌握追赶法,能利用追赶法求解线性方程组。③理解向量和矩阵范数定义,性质并掌握其计算方法.编程要求:利用Gauss消元法,追赶法解线性方程组。分析误差。 计算算法:①Gauss消元法:1. 消元过程设a?0,对k?1,2,?,n?1计算(k)kk(k)(k)?mik?aik/akk?(k?1)(k)(k)?aij?mikakj?aij?(k?1)(k)(k)?bi?mikbk?bii,j?k?1,k?2,?,n⒉回代过程(n)(n)?xn?bn/ann?n?(i)(i)(i)x?(b?ax)/a?iiijjii?j?i?1?i?n?1,?,2,1②追赶法:1.分解Ax=f: ??c/b,??c/(b?a?)( i?2,3,?,n?1)2.解Lg=f,求g:g?f/b,g?(f?ag)/(b?a?)
(i?2,3,?,n)3.解Ux=g,求x:x?y,x?y??x
(i?n?1,n?2,?,2,1) ③范数:111iiiii?1111iiii?1iii?1nniiii?1常用向量范数有:(令x=( x1,x2,…,xn)T)1-范数:
║x║1=│x1│+│x2│+…+│xn│2-范数:
║x║2=(│x1│+│x2│+…+│xn│)2221/2∞-范数: ║x║?=max(│x1│,│x2│,…,│xn│)常用的三种向量范数导出的矩阵范数是:1-范数:║A║1= max{║Ax║1/║x║1=1}=maxn1?j?n?i?1ajj2-范数:║A║2=max{║Ax║2/║x║2=1}=?1,λ1是ATA的最大特征值.∞-范数:║A║?=max{║Ax║?/║x║?=1}=maxn1?i?n?j?1aij实验例题⑴:用Gauss消元法解方程组?2???1?1??1?231??x1??3??x2?1???x3??4???????5? ??6????实验例题⑵: 用追赶法解三对角方程组Ax=b,其中?2??1?A??0??0??0?12?1000?12?1000?12?10??1????00???0?,b??0?????1??0??2???0??实验例题⑶:设?0.6A???0.1?0.5??, 0.3??计算A的行列范数.程序①:Gauss消元法function x=Gauss(A,b)%A是线性方程组的系数矩阵,b为自由项. n=length(A) for k=1:n-1m(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m(k+1:n,k)*A(k,k+1:n);
b(k+1:n)=b(k+1:n)-m(k+1:n,k)*b(k); endx=zeros(n,1); x(n)=b(n)/A(n,n); for k=n-1:-1:1x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k); end x=x';disp(sprintf('k
x(k)')); for i=0:ndisp(sprintf('%d
%f ',i,x(i+1))); end数值结果:x=Gauss(A,b)n =3k
2.555556程序②:追赶法function [x,y,beta]=zhuiganfa(a,b,c,f)%a,b,c是三对角阵的对角线上的元素,f是自由项. n=length(b);beta(1)=c(1)/b(1); for i=2:nbeta(i)=c(i)/(b(i)-a(i)*beta(i-1)); endy(1)=f(1)/b(1); for i=2:ny(i)=(f(i)-a(i)*y(i-1))/(b(i)-a(i)*beta(i-1)); endx(n)=y(n); for i=n-1:-1:1x(i)=y(i)-beta(i)*x(i+1); enddisp(sprintf('k
beta(k)')); for i=0:ndisp(sprintf('%d
%f ',i,x(i+1),y(i+1),beta(i+1))); end数值结果:a=[0 -1 -1 -1 -1]'; b=[2 2 2 2 2]'; c=[-1 -1 -1 -1 0]'; f=[1 0 0 0 0]';[x,y,beta]=zhuiganfa(a,b,c,f)k
0.000000程序③:
1.列范数:function
fan=lie(A) %A为已知矩阵 n=length(A) for j=1:n
for i=1:nx(j)=x(j)+abs(A(i,j));end endfan=max(x)disp(sprintf('n
x(n)')); for i=0:ndisp(sprintf(' %d
%f',i,x(i+1))); end 数值结果:
fan=lie(A)
fan =0.8000n
x(n) 0 0. 0.8000002.行范数:function fan=hang(A) %A为已知矩阵 n=length(A) for i=1:n
for j=1:nx(i)=x(i)+abs(A(i,j));
end endfan=max(x)disp(sprintf('n
x(n)')); for i=0:ndisp(sprintf(' %d
%f',i,x(i+1))); end数值结果:
fan=hang(A)
fan =1.1000
x(n) 0 1. 0.400000总结:从代数上看,直接分解法和Gauss消去法本质上一样,但如果我们采用”双精度累加”计算?aibi,那么直接三角分解法的精度要比Gauss消去法为高.求线性方程组的直接法,其算式繁杂,给人以枯燥沉闷的感觉.为了改善教学效果,本章着重介绍了三对角方程组的追赶法.三对角方程组以及其拓广形式的带状方程组有着广泛的实际应用.追赶法是解三对角线方程组(对角元占优势)的有效方法,它具有计算量少,方法简单,算法稳定等优点,具有鲜明的对称美.复习思考题1. 用消去法解线性方程组为什么最好选主元?怎样的方程组可以不用选主元?2. 用高斯—约当消去法求矩阵的逆,其理论根据是什么?3. 求矩阵A的LU分解的紧凑格式有何规律?写出LU分解法解Ax = b的算法。4. 乔累斯基分解法适用于哪类线性方程组?有何优点?
5. 追赶法有何优点?哪类方程组可用追赶法求解?6. 向量范数与矩阵范数是如何定义的?各有何性质?常用范数有哪些?7. 什么叫病态方程组?什么叫矩阵的条件数?条件数刻画了矩阵的什么性质?能否用选主元的方法求得病态方程组的高精度的解?
范文四:题目:煤层瓦斯含量规律分析算法:线性方程组的直接接法
组号:22组员:张玉柱 薛洪来 孔杰 商鹏煤层瓦斯含量规律分析张玉柱,薛洪来,孔杰, 商鹏(河南理工大学 安全学院,河南 焦作454000)摘要:通过煤层瓦斯含量预测数学模型的建立,研究对煤层瓦斯含量预测影响的煤层顶底板标高、埋藏深度、上覆岩层厚度、挥发分等因素,确定影响瓦斯含量的多元回归方程,为瓦斯含量的预测和估计提供了一定的理论依据。 关键词:瓦斯含量;模拟测试Mathematical models of gas content predictYuzhu-zhang, Honglai-xue, jie-kong, peng-shang(School of Safety Science and Engineering, Henan Polytechnic University,Jiaozuo.454000,China)Abstract: Through the gas content prediction mathematical model. Research
on the impact of coal seam gas content prediction top bottom elevation, the buried depth, thickness of the overlying strata, volatile gradation factors. Determine the impact of gas content and multiple regression equation for the gas content prediction and estimate provided theoretical basis.Key words:teetonicslly coal; simulation test 0. 问题背景瓦斯是指在煤矿生产过程中,从煤层、岩层和采空区放出的各种有害气体的总称,其中甲烷是瓦斯的主体成分,所以狭义的矿井瓦斯一般是指甲烷,主要来自煤层,它构成威胁煤矿开采的主要危险。它对矿井安全的威胁主要有突出、爆炸、和窒息三种形式,最严重的瓦斯灾害是瓦斯爆炸和瓦斯突出事故,它严重威胁着井下人员的生命和矿井设施的安全[1]。 瓦斯含量是影响煤矿安全生产的重要因素,因此,加强煤层瓦斯含量预测方法及瓦斯涌出的影响因素研究,掌握煤层瓦斯含量预测规律,对改善我国煤矿安全生产状况具有积极的意义[2]。本文收集、整理和分析了大量实测数据资料,通过实测和数学方法,研究对瓦斯含量预测影响的煤层顶底板标高、埋藏深度、上覆岩层厚度、挥发分等因素,确定影响瓦斯含量的多元回归方程。最后,运用该方法对夏店煤矿回采工作面进行了瓦斯含量预测,结果与现场实测数据基本吻合。根据夏店煤矿的生产实际,对瓦斯含量影响因素进行分析,研究瓦斯含量与煤层顶底板标高、埋藏深度、上覆岩层厚度、挥发分等因素之间的关系,对夏店煤矿防治矿井瓦斯灾害,确保煤矿安全生产具有重要意义。1. 问题分析及模型影响煤层瓦斯含量预测的有煤层顶底板标高、埋藏深度、上覆岩层厚度、挥发分等因素。如何判断这些因素对煤层瓦斯含量的影响程度,建立对煤层瓦斯含量影响程度大的因素的相关方程,是预测煤层瓦斯含量的关键。通过测量收集煤层瓦斯含量影响因素的实际值,代入相关方程就可以预测煤层瓦斯含量。瓦斯含量预测过程常用回归分析方法,用来研究煤层瓦斯含量y与其影响煤层瓦斯含量变量x之间的关系。简单的说,就是把对煤层瓦斯含量y有显著影响的自变量x逐个引入回归方程。首先选出与y相关程度最大的自变量,通过统计检验,表明该自变量的作用显著时,则将其引入回归式。然后在剩余的自变量中再挑选与y关系密切的自变量。当引入的自变量由于后来变量的引入使它对因变量的作用有显著变为不显著时,则随时将他们从回归式中剔除。因此,逐步回归每次在引入新变量的显著变量之前和剔除不显著变量之后,回归式中只包含显著变量。如此反复进行,直到再也没有一个自变量可以剔除为止。最终确定煤层瓦斯含量预测回归方程。确定回归方程后,把测量和收集的相应自变量数据代入。通过线性方程组Gauss消去法求得相应自变量系数,建立煤层瓦斯含量预测方程。[3]2. 多元线性回归方程原理多元线性回归分析方法,通常是把所有m 个自变量全部引入回归方程,不管自变量对因变量是否有显著影响。实际上,我们经常需要从许多自变量中“挑选”出有意义的自变量,建立一个“最优”的回归方程。所谓“最优”的回归方程,就是包含所有对因变量有显著作用的量,而不包含对因变量无显著影响的变量的回归方程。因此,在线性回归分析中发展了对自变量进行筛选的数学方法,而逐步回归分析中发展了对自变量进行筛选的数学方法,而逐步回归方程目前被认为是一种较合理的方法。回归分析是瓦斯含量预测过程一种常用的统计分析方法,用来研究某一变量与其他若干变量之间的关系,也可以定量确定一个变量随其变量变化如何变化。通常被回归的变量用y表示,称其因变量;影响因变量的其它因素用x1,x2,,,,xm表示,称为自变量。回归分析中,假设进行n次观测,观测值为(x1k,x2k,,,,xmk;yk)
(k=1,2,,, n)那么:n[4]lij??(xk?1nik?xi)(yjk?yj)??(i,j?1,2,?,m)(2-1)li0??(xk?1ik?x)(yk?y)??(i,j?1,2,?,m)(2-2)l00=错误!未找到引用源。
(2-3)y?1n?nyk(2-4)错误!未找到引用源。k?1xi均?1n?x
错误!未找到引用源。(i,j=1,2,,,m)ikk?1n(2-5)回归方程为:y?=b0?b1x1?b2x2???bmxm^(2-6)式中回归系数b1,b2,,,,,,bm由下方程组通过Gauss消去法解出:?l11b1?l12b2?..?l1mbm?l10??l21b1?l22b2?...?l2mbm?l20??..........................................?lb?lb?...?lb?lm22mmmm0?m11(2-7)常数项mb0= y??bk?1ixi(2-8)然后用方差分析法对回归方程进行显著性效果检验,说明因变量与自变量密切的相关程度。3. 算法的MATLAB实现3.1 最优回归方程和实验数据由影响夏店煤矿瓦斯含量的重要因素,把煤层底板标高、埋藏深度和上覆基岩厚度作为因变量,瓦斯含量作为因变量(见表1),进行多元回归分析。整理、计算已采区的瓦斯实际含量数据及对应的顶板泥岩和上覆基岩厚底,选取13组数据。选择瓦斯含量y为因变量,煤层厚度x1、埋藏深度x2挥发分x3和底板标高x4三个自变量(表3)建立多元线性回归方程(2-6),采用逐步回归的方法对自变量进行筛选直到建立最优的回归方程。y?b0?b1x1?b2x2?b3x3?b4x4
(3-1)式中:b0——常数项b1,b2,b3,b4——待定系数回归方程的显著性检验见表2。F统计量的值为6.0972,大于查表临界值3.86,回归是显著的,说明瓦斯含量与埋藏深度和上覆基岩厚度两个自变量有密切的相关关系,回归方程具有实际意义。表1 瓦斯含量影响因素关系表影响因素(x) 底板标高 埋藏深度瓦斯含量(y)与影响因素(x)关系式相关系数 0.4y?-0. y?0.8 y?0.9 y?0.5 y?0.1 y?-4. y?-0.50x?13.01表2 方差分析法进行显著性效果检验上覆基岩厚度0.2上覆30m内泥岩厚度下伏20m内泥岩厚度煤层厚度 煤的变质程度 (挥发份/%)0.52720.179 0.203Table 2 variance analysis significant effectinspection方差来源平方和m均方S回m统计量F?置信限Fa统计推断 当F>Fa时,认为回归显著,线性相关密切;当F>Fa时,认为回S回?回归?bll?1ii0s回?s回s余126.3049S余?l00?S回42.1016s余?S余n?m?1(m,n-m-1)6.09723.86归不显著,线性相关不密切。其中Fa=3.86剩余39.62434.4027表3 夏店煤矿3号煤层瓦斯含量及其影响因素统计表Table 3
No. 3 coal mine of Xiadian factors affecting gas content statistical table上覆30m内泥岩厚度下伏20m内泥岩厚度 7.6挥发份上覆基岩厚度样品 煤厚 埋深 顶板岩性 标高 瓦斯含量JJ-1 6.1 435.8 15.3 8.55 砂质泥岩 421.26 491.748 5.5ZK3-1 ZK3-35.9 5.65427.6 523.814.85 12.8 上覆30m内泥岩厚度12.55 15 下伏20m内泥岩厚度 6.29 6.5 10.33 9.5511.41 9.62 挥发份泥岩 砂质泥岩400.8 517.55 上覆基岩厚度565.261 408.6353.2 13样品 煤厚 埋深 顶板岩性 标高 瓦斯含量ZK1-1 ZK1-2 ZK6-2 ZK7-1 503 604 602 804 802 7095.54 5.6 5.6 6.2 5.9 6.32 5.72 6.15 5.23 5.93290.67 368.75 426.8 386.95 188.41 433.81 282.84 536.47 288.93 480.757.98
17.37 14.9510.44 9.78 9.98 15.42 11.35 13.29 10.27 9.73 12.02 11.87砂质泥岩 砂质泥岩 泥岩 泥岩 泥岩 泥岩 泥岩 泥岩 泥岩 砂质泥岩268.3 358.65 417.4 360.4 147.02 420.24 274.33 483.83 271.01 471.25639.952 588.895 500.607 596.527 714.98 458.87 675.56 458.71 640.15 441.231.7 1.3 7 5.2 1 6.7 3.3 11.3 4.1 8.33.2 Matlab程序代码确定估计系数b0, b1, b2, b3,b4原则是根据线性方程组Gauss消去法[5],把表3数据输入Matlap程序如下:function b=gauss(A,x) % 高斯求解方程组 %b=gauss(A,x) n=length(A); a=[A,x]; for k=1:n-1maxa=max(abs(a(k:n,k))); if maxa==0 endfor i=k:nif abs(a(i,k))==maxay=a(i,k:n+1);a(i,k:n+1)=a(k,k:n+1);a(k,k:n+1)=y; end endfor i=k+1:nl(i,k)=a(i,k)/a(k,k);a(i,k+1:n+1)=a(i,k+1:n+1)-l(i,k).*a(k,k+1:n+1); endend%回代if a(n,n)==0 return endx(n)=a(n,n+1)/a(n,n); for i=n-1:-1:1x(i)=(a(i,n+1)-sum(a(i,i+1:n).*x(i+1:n)))/a(i,i); end调用示例如下:>>A=[1,6.1,435.8,8.55,491.748;1,5.9,427.6,11.41,565.261;1,5.65,523.8,9.62,408.635;1,5.54,290.67,10.44,639.952;1,6.2,386.95,15.42,596.527;1,5.72,282.84,10.27,675.56;1,6.15,536.47,9.73,458.71;1,5.93,480.75,11.87,441.23]; >> x=[5.5; 3.2; 13; 1.7; 5.2;3.3;11.3;8.3]; >>b=gauss(A,x)b=18.5 0.7 -0.0209可求得:b0=18.2583,b1=-1.3535,b2=0.1350,b3=0.1267,b4=-0.0209 则回归方程为:y=18.5x1+0...0209x4
(3-2)4. 煤层瓦斯含量预测结果及分析4.1 煤层瓦斯含量预测结果通过以上分析,得出该矿影响瓦斯含量的多元回归方程如下,预测井田北部瓦斯含量。y?18.5x1?0...0209x4(4-1)依据北部瓦斯资料,统计了煤层顶底板标高(x)与瓦斯含量(y)(见表4)并进行线性回归分析,煤层底板标高一定程度上反映了煤层埋厚的变化,分析线性相关性较好,故依此预测井田北部煤层瓦斯含量:结果如下:表4 夏店煤矿3#煤层瓦斯含量预测结果表Table 4
No. 3 coal mine of Xiadian gas content prediction results table瓦斯含量趋势值 南部底板标高/m2.1 +6823.2 +6494.1 +6175.6 +5846.5 +5527.4 +5208.3 +4879.1 +45510.2 +44211.8 +3994.2 煤层瓦斯含量预测结果分析井田北部瓦斯含量实测值(见表5)与煤层瓦斯含量预测值进行比较,误差均在允许范围内。说明这种瓦斯含量预测方法正确合理。表5 夏店煤矿3#煤层瓦斯含量实测值Table 3
No. 3 coal mine of Xiadian gas content actual measured value瓦斯含量趋势值 南部底板标高/m2.1 +6823.1 +6494.2 +6175.6 +5846.6 +5527.3 +5208.3 +4879 +45510.1 +44211.8 +3995. 结论通过对煤层瓦斯含量进行多元线性回归分析,找出了煤层瓦斯含量随煤层底板标高、埋藏深度和上覆基岩厚度变化规律。利用线性方程gauss消去法对回归方程进行求解,确定了煤层瓦斯含量与煤层底板标高、埋藏深度和上覆基岩厚度变化方程本文所测定的瓦斯含量规律是夏店煤矿3#工作面瓦斯涌出的一般规律,对不同煤层工作面瓦斯含量规律,需要测定有关参数来修正。参考文献[1] 林柏泉,崔恒信.矿井瓦斯防治理论与技术[M].徐州:中国矿业大学出版社,1998. [2] 林井祥,孙广义.煤与瓦斯突出的突出规律的参数研究[J].矿冶—2011年3期.[3] 魏风清,张建国.煤与瓦斯突出预测方法和防止措施[M].北京:煤炭工业出版社,2005. [4] (美)韦斯伯格,王静龙.应用线性回归[M].北京:中国统计出版社,1998.3. [5] 李庆扬,王能超,易大义.数值分析[M].武汉:华中科技大学出版社,2006.
范文五:数值分析第一次上机实习报告——线性方程组直接解法一、问题描述设 Hn = [hij ] ∈ Rn×n
是 Hilbert 矩阵, 即hij?1 i?j?1对n = 2,3,4,…13,?1???(a) 取x????Rn?n,及bn?Hnx,用Gauss消去法和Cholesky分解方法来求解?1???Hny?bn,看看误差有多大.(b) 计算条件数:cond(Hn)2(c) 使用某种正则化方法改善(a)中的结果.二、方法描述1. Gauss消去法Gauss消去法一般用于系数矩阵稠密且没有任何特殊结构的线性方程组。设H=[hij],y = (y1,y2,…,yn)T. 首先对系数矩阵Hn进行LU分解,对于k=1,2,…n,交替进行计算:k?1??ukj?hkj??lkrurj),j?k,k?1,…,n?r?1 ?k?11?l?(aik??lirurk),i?k?1,k?2,…,nik?ukkr?1?由此可以得到下三角矩阵L=[lij]和上三角矩阵U=[uij]. 依次求解方程组Ly=b和Ux=y ,i?1??yi?bi??liryr,i?1,2,…,n?r?1 ?n1?x?(y?uirxr),i?n,n?1,…,1?ii?uiir?i?1?即可确定最终解。2. Cholesky分解法对于系数矩阵对称正定的方程组Hny?bn,存在唯一的对角元素为正数的下三角矩阵L,使得H=LLT。因此,首先对矩阵Hn进行Cholesky分解,即1j?1?22?ljj?(hjj??ljk)k?1?
i?j?1,…n ?j?11?l?(h?ll)?ijijikjk?lk?1jj?L的元素求出之后,依次求解方程组Ly=b和LTx=y,即b1?y??1l?11 ?i?11?y?(b?ly),i?2,3,…n?iiikk?lk?1ii?yn?x??nl?nn ?n?x?1(y?lkixk),i?n?1,n?2,…1?ii?lnnk?i?1?由此求得方程组Hny?bn的解。3. 病态方程组预处理通过对病态方程组进行预处理,可以降低系数矩阵的条件数。选择Hn的对角元素进行开方构成对角矩阵D,从而构造新的系数矩阵G,使得G= D-1HnD-1。从而使得原方程Hny?bn,转化为等价方程组Gz?D-1b则原方程的解y?D?1z。三、方案设计通过MATLAB进行编程,对上述问题进行求解。在所有的程序文件中,main.m为主程序文件,即对方程组采用Gauss消去法和Cholesky分解法分别进行求解,并计算不同阶数Hilbert 矩阵的条件数cond(Hn)2;gauss.m文件为Gauss消去法解方程的matlab算法,其中调用了两个函数,分别为求解下三角形方程组的ltri函数和求解上三角形方程组的utri函数,这两个函数相应的分别编写在ltri.m文件和utri.m文件中。通过平衡法对Hilbert 矩阵进行预处理,调用文件balance.m中的算法进行计算和处理。四、计算结果及其分析a) 用Gauss消去法和Cholesky分解方法来求解方程Hny?bn,结果如下b) 计算条件数cond(Hn)2,结果如下c) 通过对矩阵进行预处理,改善a) 中的计算结果。
范文六:第4章线性方程组的数值解法直接法迭代法2.求解方法(1)直接解法:经过有限次的算术运算(四则运算),计算出方程组的精确解(在每一步运算过程都不会发生舍入误差的假设下)。但实际计算中不可避免舍入误差的存在和影响,所以只能求得线性方程组的近似解。适用于低阶(n?1000)稠密(零元素较少)的方程组及系数矩阵为带状的方程组。主要方法有:1).Gramer法则:实际计算中不常用,可用于理论分析。优点:理论上得到的解是精确的缺点:计算量大,共(n?1)(n?1)n!?n次乘、除运算,占用空间大2).Gauss消去法:逐次消元,再回代求解,实际计算中常用。优点:计算量相对Gramer方法小很多,便于程序化缺点:误差分析困难,当用很小的数作除数或两个相近数相减时,往往造成其解不稳定3).改进的Gauss消去法:实际计算中常用的有效方法。如:列(行)主元法,全主元法,LU分解法,平方根法,追赶法等(2)迭代法:利用不动点理论将方程组变形为某种等价的迭代公式,给出初始解x,求出迭代序列x(k?0,1,2,?)的极限x(精确解)。适宜于求解大型(n?1000)稀疏的(零元素较多)线性方程组。*(0)(k)主要方法有:1).雅可比(Jacobi)迭代法2).高斯—赛德尔(Guass?Seidel)法是G?S方法的加速3).超松弛(SOR)法:优点:速度高,计算机存储量小,程序设计简单;缺点:有误差,需要事先确定迭代算法的收敛性及收敛速度。(3)算法的评价:评价求解Ax?b数值方法好坏的标准:1).近似解与真解的误差大小(常用范数度量)2).运算量或计算时间(加、减、乘、除的运算次数或步数,通常只记乘、除法)3).存储空间4).逻辑复杂度§1
Gauss消去法1. Gauss消去法: 设有线性代数方程组?a11?a?21????an1a12?a1n??x1??b1??????a22?a2n??x2??b2????????????
(1.1)?????an2?ann??xn??bn?A?(aij)n?n,b?(bi)n?1记作
Ax?b,基本思想: 首先通过逐次消元将Ax?b化为等价的上三角形方程组(称为消元过程);再求解上三角方程组(称为回代过程)。?0.x2?2例如方程组?x1?x2?3?准确到小数点后第9位的解为x1?2.,x2?0..如果计算过程用四位十进制浮点数(仿机器实际计算),用第一个方程消去第二个方程中的,得10?0.?0.?0.2000??66?10?0.?0.2000?由此解得x2?1,x1?0,显然它不是原方程的解.?4为了减少计算过程中的舍入误差对解的影响,应选择绝对值尽可能大的主元作除数,基于这种思想导出了选主元消去法。?2x1?x2?2x3?5例:求解线性方程组??5x1?x2?x3?8??x1?3x2?4x3??4?212??5??5?11?8?x3?2,x2??1,???4??1?3?4???第1列的列主元a21?5,回代求解?2行与1行交换?5?11?8????212?5??5?11??0?2.8?4.2??1?3?4??4??????00?0.5??第1行分别乘以?2/5,?1/5并加到2,3行2行乘以1/2加到3行???5?11?8???01.41.6?1.8?????????第2列的列主元a32??2.8,?5?11??2.8?4.2??5.6??0?2.8?4.2??0???3行与2行交换??01.41.6?x1?18??5.6??1???8??5.6?1.8???注意:列交换改变了xi的顺序,需记录交换的次序,解完后再交换回来。为方便起见,计算时可以做一个存储表。4. Gauss-Jordan消去法(无回代的Gauss消去法) 基本思想: 在第k步选取主元a(k)kk(k)kk?0后,将第k列其它元全消为零元,再将a化为1。同解方程为(1)x?b1??10...0??1??1...0??x??(2)?b2??????2??????????? ?????(n)?1??xn??bn???2x1?x2?2x3?5?例:用高斯?若当消去法求解线性方程组?5x1?x2?x3?8?x?2x?4x??423?1?212??5?11??1?3?4??5?11??212??1?3?4???5?11?01.41.6??0?2.8?4.25?8??4???8?5??4???8?1.8??5.6????5005???1?0?2.802.8????0??00?0.5?1?????0???502.510??0?2.8?4.2?5.6???0.5?1??00?????5?118???0?2.8?4.2?5.6??1.4?0.51.8??0??001?10?1?012????x1?1,x2??1,x3?2高斯?若当法计算复杂度:解方程组大约需要n/2次乘除法,比高斯消去法大,因此一般并不用于解方程组,而是矩阵求逆.高斯?若当法求矩阵的逆:对分块矩阵?A?I?作高斯?若当消元,3当把A化为单位阵I时,单位矩阵I即化为A?1§2
矩阵的LU分解及对称正定矩阵的平方根法下面我们进一步利用矩阵理论来分析Gauss消去法。1.Gauss消去法的矩阵形式已知A?R,并假设在Gauss消元法中an?n(k)kk?0(k?1:n),(1)(1)(1)(2)(2)(2)B?(A,b)?(A,b)?B则将 ,相当于?1?(1)(1)???a11a12??l211??a(1)a(1)??l31??21221???????????a(1)a(1)??l??n1n21?n1?(1)(1)?a?a??(1)1n(1)2n(1)?ann(1)(1)?b?a11a12??0a(2)b??22???????(1)?(2)?bn???0an2(1)1(1)2?a?a???a(1)1n(2)2n(2)nn????
(2)?bn??bb?(1)1(2)2记作
L1(A,b)?(A,b)(2)(2)(k)(k)(k)(k?1)(k?1)(k?1)B?(A,b)?(A,b)?B将,相当于?1????1??lk?1,k1??????ln,k?(1)(1)(1)?a11a?a1?1,k?1k???????(k)(k)?(k)?akkak,k?1?B??(k?1)?a?k?1,k?1??????(k?1)?1?a?n,k?1???????(1)a1n?(k)ak,n(k?1)ak?1,n?(k?1)an,nb1(1)????????B(k?1)(k?1)?bk?1???(k?1)?bn?(k)(k)(k?1)(k?1)L(A,b)?(A,b) 记作
k经过n-1(1)(1)(n)(n)LL?LL(A,b)?(A,b),即步后得
n?1n?221(n)???LLLLAA?n?1n?221?(n) ??Ln?1Ln?2?L2L1b?b记L?LL?L?11?12?1n?1?1?l?21??l31????l?n11l32?ln21?ln3??u11????,U?A(n)?????????1??u12u22u1n??...u2n?????unn?...于是A?LU——称为矩阵A的LU(三角)分解。显然L是单位下三角矩阵,U是上三角矩阵。由于A?LU,则求解方程组Ax?b?求解下列两个三?Ly?b角形线性方程组? ?Uxy?计算公式分别为?y1?b1?i?1?及???yblb(i2:n)?iiijj?j?1??xn?yn/unn?n?x?(y?ux)/u(i?n?1:?1:1) ?iiijjii?j?i?1?问题: 一般的,矩阵A不一定能够作LU分解。那么,在什么条件下能保证A可作LU分解?并且唯一?例
用Gauss?x1?x2?x3?6?法的LU分解求解方程组?4x2?x3?5?2x?2x?x?123?1解 由Gauss消元法 l21?0,l31?2,l32??1,?100??111?????
A??010??04?1??LU??2?11????00?2??求解Ly?b,得y1?6,y2?5,y3??6, 求解Ux?y,得x3?3,x2?2,x1?1。2. 矩阵的直接LU分解法实际上也前面由Gauss消去法给出了A的LU分解。可以直接从矩阵A的元素得到计算L,U元素的递推公式,而不需要任何中间步骤,这就是所谓的矩阵A的LU直接分解法。由表可知??1??1A??01??????121???0101????解Ly?b,得y?(5,3,6,4)T解Ux?y,得x?(1,1,2,2)T。 02020?1?1??2??1定理2
若A?R是对称矩阵,且A的顺序主子式Di?0(i?1:n),则A可唯一分解为n?nA?LDL。其中L为单位下三角矩阵,D为对角矩阵。证明:因为A的各阶顺序主子式都不等于零,由定理1知,A可惟一分解为T?1?l21?A?????ln11?ln2??u11??????????1??u12u22?u1n???u2n?且uii?0(i?1,2,?,n)????unn?3.
实对称正定矩阵的平方根法应用有限元法解结构力学的问题时,最后归结为求一个线性方程组,其系数矩阵大多是对称正定矩阵。下面我们研究对称正定矩阵的分解方法定义
设A为n阶实对称矩阵,若对任意的非零向量x,都有xAx?0 则称A为对称正定矩阵。T定理3
设A为对称正定矩阵,则下列命题等价 (1)A是正定矩阵;(2)A可逆,且A亦为对称正定;?1(3)A的任意主子阵Ai(i?1:n)亦为对称正定; (4)A的特征值?i?0;(5)A的顺序主子式Di?0(i?1:n)。(1) Cholesky—乔列斯基分解 定理
设A?Rn?n对称正定矩阵,则存在唯一的主对角线元素都是正数的下三角矩阵L,使得A?LL?a11a12?aa2122??即 ????an1an2?a1n??l11???a2n??l21?????????ann??ln1l22?ln2??l11??????????lnn??l21?ln1??l22?ln2?????lnn?T称该式为对称正定矩阵A的Cholesky分解。(2)平方根法若将对称正定矩阵A分解为A?LL,则求解方程组TAx?b?求解下列两个三角形线性方程组?Ly?b?TLx?y?称此方法为平方根法。第4章线性方程组的数值解法直接法迭代法2.求解方法(1)直接解法:经过有限次的算术运算(四则运算),计算出方程组的精确解(在每一步运算过程都不会发生舍入误差的假设下)。但实际计算中不可避免舍入误差的存在和影响,所以只能求得线性方程组的近似解。适用于低阶(n?1000)稠密(零元素较少)的方程组及系数矩阵为带状的方程组。主要方法有:1).Gramer法则:实际计算中不常用,可用于理论分析。优点:理论上得到的解是精确的缺点:计算量大,共(n?1)(n?1)n!?n次乘、除运算,占用空间大2).Gauss消去法:逐次消元,再回代求解,实际计算中常用。优点:计算量相对Gramer方法小很多,便于程序化缺点:误差分析困难,当用很小的数作除数或两个相近数相减时,往往造成其解不稳定3).改进的Gauss消去法:实际计算中常用的有效方法。如:列(行)主元法,全主元法,LU分解法,平方根法,追赶法等(2)迭代法:利用不动点理论将方程组变形为某种等价的迭代公式,给出初始解x,求出迭代序列x(k?0,1,2,?)的极限x(精确解)。适宜于求解大型(n?1000)稀疏的(零元素较多)线性方程组。*(0)(k)主要方法有:1).雅可比(Jacobi)迭代法2).高斯—赛德尔(Guass?Seidel)法是G?S方法的加速3).超松弛(SOR)法:优点:速度高,计算机存储量小,程序设计简单;缺点:有误差,需要事先确定迭代算法的收敛性及收敛速度。(3)算法的评价:评价求解Ax?b数值方法好坏的标准:1).近似解与真解的误差大小(常用范数度量)2).运算量或计算时间(加、减、乘、除的运算次数或步数,通常只记乘、除法)3).存储空间4).逻辑复杂度§1
Gauss消去法1. Gauss消去法: 设有线性代数方程组?a11?a?21????an1a12?a1n??x1??b1??????a22?a2n??x2??b2????????????
(1.1)?????an2?ann??xn??bn?A?(aij)n?n,b?(bi)n?1记作
Ax?b,基本思想: 首先通过逐次消元将Ax?b化为等价的上三角形方程组(称为消元过程);再求解上三角方程组(称为回代过程)。?0.x2?2例如方程组?x1?x2?3?准确到小数点后第9位的解为x1?2.,x2?0..如果计算过程用四位十进制浮点数(仿机器实际计算),用第一个方程消去第二个方程中的,得10?0.?0.?0.2000??66?10?0.?0.2000?由此解得x2?1,x1?0,显然它不是原方程的解.?4为了减少计算过程中的舍入误差对解的影响,应选择绝对值尽可能大的主元作除数,基于这种思想导出了选主元消去法。?2x1?x2?2x3?5例:求解线性方程组??5x1?x2?x3?8??x1?3x2?4x3??4?212??5??5?11?8?x3?2,x2??1,???4??1?3?4???第1列的列主元a21?5,回代求解?2行与1行交换?5?11?8????212?5??5?11??0?2.8?4.2??1?3?4??4??????00?0.5??第1行分别乘以?2/5,?1/5并加到2,3行2行乘以1/2加到3行???5?11?8???01.41.6?1.8?????????第2列的列主元a32??2.8,?5?11??2.8?4.2??5.6??0?2.8?4.2??0???3行与2行交换??01.41.6?x1?18??5.6??1???8??5.6?1.8???注意:列交换改变了xi的顺序,需记录交换的次序,解完后再交换回来。为方便起见,计算时可以做一个存储表。4. Gauss-Jordan消去法(无回代的Gauss消去法) 基本思想: 在第k步选取主元a(k)kk(k)kk?0后,将第k列其它元全消为零元,再将a化为1。同解方程为(1)x?b1??10...0??1??1...0??x??(2)?b2??????2??????????? ?????(n)?1??xn??bn???2x1?x2?2x3?5?例:用高斯?若当消去法求解线性方程组?5x1?x2?x3?8?x?2x?4x??423?1?212??5?11??1?3?4??5?11??212??1?3?4???5?11?01.41.6??0?2.8?4.25?8??4???8?5??4???8?1.8??5.6????5005???1?0?2.802.8????0??00?0.5?1?????0???502.510??0?2.8?4.2?5.6???0.5?1??00?????5?118???0?2.8?4.2?5.6??1.4?0.51.8??0??001?10?1?012????x1?1,x2??1,x3?2高斯?若当法计算复杂度:解方程组大约需要n/2次乘除法,比高斯消去法大,因此一般并不用于解方程组,而是矩阵求逆.高斯?若当法求矩阵的逆:对分块矩阵?A?I?作高斯?若当消元,3当把A化为单位阵I时,单位矩阵I即化为A?1§2
矩阵的LU分解及对称正定矩阵的平方根法下面我们进一步利用矩阵理论来分析Gauss消去法。1.Gauss消去法的矩阵形式已知A?R,并假设在Gauss消元法中an?n(k)kk?0(k?1:n),(1)(1)(1)(2)(2)(2)B?(A,b)?(A,b)?B则将 ,相当于?1?(1)(1)???a11a12??l211??a(1)a(1)??l31??21221???????????a(1)a(1)??l??n1n21?n1?(1)(1)?a?a??(1)1n(1)2n(1)?ann(1)(1)?b?a11a12??0a(2)b??22???????(1)?(2)?bn???0an2(1)1(1)2?a?a???a(1)1n(2)2n(2)nn????
(2)?bn??bb?(1)1(2)2记作
L1(A,b)?(A,b)(2)(2)(k)(k)(k)(k?1)(k?1)(k?1)B?(A,b)?(A,b)?B将,相当于?1????1??lk?1,k1??????ln,k?(1)(1)(1)?a11a?a1?1,k?1k???????(k)(k)?(k)?akkak,k?1?B??(k?1)?a?k?1,k?1??????(k?1)?1?a?n,k?1???????(1)a1n?(k)ak,n(k?1)ak?1,n?(k?1)an,nb1(1)????????B(k?1)(k?1)?bk?1???(k?1)?bn?(k)(k)(k?1)(k?1)L(A,b)?(A,b) 记作
k经过n-1(1)(1)(n)(n)LL?LL(A,b)?(A,b),即步后得
n?1n?221(n)???LLLLAA?n?1n?221?(n) ??Ln?1Ln?2?L2L1b?b记L?LL?L?11?12?1n?1?1?l?21??l31????l?n11l32?ln21?ln3??u11????,U?A(n)?????????1??u12u22u1n??...u2n?????unn?...于是A?LU——称为矩阵A的LU(三角)分解。显然L是单位下三角矩阵,U是上三角矩阵。由于A?LU,则求解方程组Ax?b?求解下列两个三?Ly?b角形线性方程组? ?Uxy?计算公式分别为?y1?b1?i?1?及???yblb(i2:n)?iiijj?j?1??xn?yn/unn?n?x?(y?ux)/u(i?n?1:?1:1) ?iiijjii?j?i?1?问题: 一般的,矩阵A不一定能够作LU分解。那么,在什么条件下能保证A可作LU分解?并且唯一?例
用Gauss?x1?x2?x3?6?法的LU分解求解方程组?4x2?x3?5?2x?2x?x?123?1解 由Gauss消元法 l21?0,l31?2,l32??1,?100??111?????
A??010??04?1??LU??2?11????00?2??求解Ly?b,得y1?6,y2?5,y3??6, 求解Ux?y,得x3?3,x2?2,x1?1。2. 矩阵的直接LU分解法实际上也前面由Gauss消去法给出了A的LU分解。可以直接从矩阵A的元素得到计算L,U元素的递推公式,而不需要任何中间步骤,这就是所谓的矩阵A的LU直接分解法。由表可知??1??1A??01??????121???0101????解Ly?b,得y?(5,3,6,4)T解Ux?y,得x?(1,1,2,2)T。 02020?1?1??2??1定理2
若A?R是对称矩阵,且A的顺序主子式Di?0(i?1:n),则A可唯一分解为n?nA?LDL。其中L为单位下三角矩阵,D为对角矩阵。证明:因为A的各阶顺序主子式都不等于零,由定理1知,A可惟一分解为T?1?l21?A?????ln11?ln2??u11??????????1??u12u22?u1n???u2n?且uii?0(i?1,2,?,n)????unn?3.
实对称正定矩阵的平方根法应用有限元法解结构力学的问题时,最后归结为求一个线性方程组,其系数矩阵大多是对称正定矩阵。下面我们研究对称正定矩阵的分解方法定义
设A为n阶实对称矩阵,若对任意的非零向量x,都有xAx?0 则称A为对称正定矩阵。T定理3
设A为对称正定矩阵,则下列命题等价 (1)A是正定矩阵;(2)A可逆,且A亦为对称正定;?1(3)A的任意主子阵Ai(i?1:n)亦为对称正定; (4)A的特征值?i?0;(5)A的顺序主子式Di?0(i?1:n)。(1) Cholesky—乔列斯基分解 定理
设A?Rn?n对称正定矩阵,则存在唯一的主对角线元素都是正数的下三角矩阵L,使得A?LL?a11a12?aa2122??即 ????an1an2?a1n??l11???a2n??l21?????????ann??ln1l22?ln2??l11??????????lnn??l21?ln1??l22?ln2?????lnn?T称该式为对称正定矩阵A的Cholesky分解。(2)平方根法若将对称正定矩阵A分解为A?LL,则求解方程组TAx?b?求解下列两个三角形线性方程组?Ly?b?TLx?y?称此方法为平方根法。
范文七:线性方程组的解法1 引言在科学研究和大型工程设计中出现了越来越多的数学问题,而这些问题往往需要求数值解。在进行数值求解时,经离散后,常常归结为求解形如Ax= b的大型线性方程组。而如插值公式,拟合公式等的建立,微分方程差分格式的构造等,均可归结为求解线性方程组的问题.在工程技术的科学计算中,线性方程组的求解也是最基本的工作之一.因此,线性 方程组的解法一直是科学和工程计算中研究最为普遍的问题,它在数值分析中占有极其重要的地位。20世纪50年代至70年 代,由于电子计算机的发展,人们开始考虑和研究在计算机上用迭代法求线性方程组Ax =b的近似解,用某种极限过程去逐渐逼近精确解,并发展了许多非常有效的迭代方法,迭代法具有需要计算机存储单元少、程序设计简单、原始系数矩阵在计算过程中始终不变等优点。例如Jacobi方法、Gauss—Seidel 方法、SOR方法、SSOR方法 ,这几种迭代方法是最常用的一阶线性定常迭代法。2 主要算法20世纪50年代至70年代,人们开始考虑和研究用迭代法求解线性方程组。Ax = b
(1)的近似解,发展了许多有效的方法,其中有Jacobi方法、Gauss—Seidel方法,SOR方法、SSOR方法,这几种迭代方法均属一阶线性定常迭代法,即若系数矩阵A的一个分裂:A =M-N ;M 为可逆矩阵,线性方程组(1)化为 :(M-N)X =b;→M X = NX + b;→X= MNX+ Mb得到迭代方法的一般公式:X(k+1)=HX(k)+d
(2)其中:H =MN-1,d=M-1b,对任意初始向量X(0)一阶定常迭代法收敛的充分必要条件是: 迭代矩H的谱半径小于1,即ρ(H)2.1 Jacobi迭代法若D为A的对角素构成的对角矩阵,且对角线元素全不为零。系数矩阵A的一个分解:A =D-(L +U); 这里D为A的对角素构成的对 角矩阵, 为严格下三角阵,U为严格上三角阵。 Jacobi迭代的矩阵形式为 :X(k+1)=HJX(k)+dJ
(3)(3) 式 中: dJ= D-1b;
HJ=I-D-1A, 称HJ为Jacobi迭代矩阵. 其计算公式为:迭代矩阵HJ的谱半径ρ(H)2.2 Gauss-Seidei迭代法对于非奇异方程组,若D为A的对角素构成的对角矩阵,且对角线元素全不为零;系数矩阵A的一个分解:A = (D-L)-U
(5)Gauss-Seide迭代矩阵形式为:
X(k+1)=HGSX(k)+dGS上式中: HGS=(D-L)-1U;dGS=(D-L)-1L称HGS为G-S迭代矩阵。计算公式为:i=1,2,3,,,n
k=0,1,2,,,
(6)若A为严格或不可约对角占优矩阵,或A为对称正定阵,则对于任意初值X(0),Gauss-Seide迭代法收敛。2.3
SOR(successive over relaxation)迭代法对于非奇异方程组,若D为A的对角素构成的 对角矩阵,且对角线元素全不为零;系数矩阵A的一个分解:?1??1?w?A??D?L???D?U? ?w??w?
(7)这里D为A的对角素构成的对角矩阵,为严格下三角阵,为严格上三角阵。SOR迭代法的矩阵形式为
X(k+1)=HSOR X(k)+dSOR式中 :HSOR=(D-wL)-1 [ (1-w) D + wU]d SOR= w ( D-wL )-1b。计算公式为:i=1,2,3,,,n
k=0,1,2,,,
(8)然后按相反的次序(i= n,n-1,,,1)用向后的SOR方法逐点计算i= n,n-1,,,1;k=0,1,,,
(11)若AX= b中A是正定的,且0﹤w﹤2,则SSOR法收敛。2.4
消元法消元法是初等数学中求解低阶多元线性方程组的方法.此时线性方组必须是适定方程组.一般是用于二元一次 或三元一次方程组.当未知元增多时.计算效率低甚至无法求解.2.5
克拉默法则当系数行列式不为零时.适定方程组有惟一 解 .其解如(12)式所示 :xi= Di /D
i=1,2,,,,n
(12)其中D是系数行列式,Di是在系数行列式基础之上结合方程组右边常数形成的新行列式。在此法则中。行列式的计算显得非常重要。用行列式的性质计算行列式最为有效。 对于二 、三阶行列式可以利用对角线法则计算。克拉默法则克服了消元法计算效率低甚至无法计算多元一次方程组的缺点.但是对于系数行列式等于零以及欠定或者超定方程组的情况,它是无能为力的。事实上,当未知元数过多时 (如未知元数≥5)。克拉默法则的计算效率就很低。2.6 逆阵乘积法对于适定方程组.可以把它表达成矩阵方程的形式:AX=b解矩阵可以利用逆阵乘积法求出:A-1 b = X矩阵运算的实质是把矩阵当作一个“量 ”来运算。使 普通数 的运算有很大的简化。但是该方法的前提是A可逆。本质上仍然是系数行列式|A|≠0.对于阶数比较高的系数矩阵.直接求解逆阵也是比较 困难的(利用初等变换可以降低求解难度)。当|A|=0时,或者对于欠定或者超定方程组,逆阵乘积法仍然是无能为力的或不适合的。2.6
初等变换法对于欠定或超定方程组的求解。初等变换法是最直接、 最简单的方法。同时该方法也能用于适定方程组的求解。因此,初等变换法是一种求解线性方程组的普适方法和最 基本方法 。秩是矩 阵的本征参 数.利用系数矩阵的秩和增广矩阵秩的关系,可以判断线性方程组解的情况:无解,惟一解和无穷多解。对增广矩阵进行初等变换.可以得到增广矩阵的等价矩阵.从而得到了原来方程组的等价方程组。由于等价矩阵已变换成阶梯形矩阵或最简形矩阵,所以等价方程组变 得非常简单。可以方便地求出解。初等变换也是求逆阵的一种直观方法.所以可以和阵乘积法配合运用。2.7
利用向量空间概念求解线性方程组如果把齐次线性方程组的解集看作是一个向量空间 。由于一个向量组 (向量空间)与它的一个最大无关组等价。则只需求出解集的一基础解系即可确定齐次线性方程组的解。求解的方法仍然采用初等变换法.解的形式可以用通解来表达.更能说明解的本质。尤其是无穷多解。基础解系中线性无关的向量形成解向量空间的一个最大无关组 。在非齐次方程组求解中.只要求出对应齐次方程组的通解。然后加上非齐次方程组的一个特解即可.而且这种通解和特解可以在一次初等变换中求出。2.7
Krylov子空间方法Krylov子空间方法是解决大型线性方程组的一类十分有效的方法 ,其主要代表是对称正定线性方程组的共轭梯度法和非对称线性方程组的GMRES方法。20世 纪 50年代 初 期 由 Hestenes和 Stiefel首先提出的共轭梯度法 (或简称 CG法 )。近20年来有关的研究取得了前所未有的发展,目前有关的方法和理论已经相当成熟,并且成为求解大型稀疏线性方程组十分有效的一类方法。通过对经典迭代法的总结 ,发现 SOR迭代法在松弛因子w取最优松弛因子wopt的情况下,迭代收敛很快。但是,计算wopt还需要求得对应的 Jacobi迭代矩阵的谱半径,这常常是比较困难的。考虑线性方程组:
(13)的求解问题,其中A是给定的n阶对称正定阵,b是给定的n维向量, x是待求的n维向量。为此定义二次泛函数Ψ( x) = xTA x -2bTx
(14)则可以证明求方程组 (13)的解 等价于求二次泛函(14)的极小点。由此,给定了初始向量x(0),按某一方向去求式(14)的极小值点x(1),就得到下一个迭代值x(1),再由x(1)出发 ,求x(2)等等。这样来逼近精确解x*.若取求最小值的方向为Ψ在x(k-1)(k =1, 2,?)处的负梯度方向,就是所谓的最速降法。 然而理论和实际计算表明这个方法的收敛速度较慢,共轭梯度法则是在x(k-1)处的梯度方向r(k-1)和
这一步的修正方向p(k-1)所构成的二维平面内,寻找使Ψ减少最快的方向作为下一步的修正方向,即求极小值的方向p(k) (其第一步仍取负梯度方向).计算公式为p(1)=r(0)=b- x(0)再逐次计算qk=( r(k-1),r(k-1))/(A p(k), p(k)),((x,y)表示向量xy的内积x(k)= x(k-1)+ qk p(k)r(k)= r(k-1)- qk p(k)ek=(r(k),r(k))/(r(k-1),r(k-1))p(k+1)= r(k)+ ek p(k)(k=1,2,,,)可以证明当i≠j时,(p(i),A p(j))=0,(ri,rj)=0从而p(1),p(2),,.形成一共轭向量组;r(0),r(1),,形成一正交向量组。后者说明若没有舍入误差的 ,至多n次迭代就可以得到线性方 程组(14)的精确解。然而,在实际计算中一般都有舍入误差,所以r(0),r(1),,并不是真正相交,而x(0),x(1),,等也只是逐步逼近式(14)的真解,故我们将共轭梯度法作为迭代法来使用。3 数值实验下面通过一个数值实例来比较Jacobi方法、Gauss-Seidel方法、SOR方法的收敛速度。 数值实验:我们取一个系数矩阵为64阶的线性方程组:A X= b。其中:精确解:X=(1
1)T1×64; X(0) (0
1)T1×64; 精度为:0.0001. 执行C语言程序interrator-method. C得到result.txt文件。得到三种迭代法的对比表1.表1 三种方法的迭代性能MethodJacobiGauss-SeidelSORwithw=1.0SORwithw=woptSORwithw=1.5 IT
CPU Time/ms 18
2.00000 ‖x(k)-x∞‖∞ Info 8.1 几乎比 Jacobi收1.68O24×10-5 敛速度快一倍 1.6 退化成 Gauss
_ Seidel迭代 6
1.552×10-6 使 SOR迭 代
达到最快 17
3.836×10-5 w的取值于SOR迭代的收敛速度影响很大4 优点克拉默法则、逆阵乘积法只能求解系数行列式不为零的适定方程组 ;初等变换法可以直观地解决所有类型的超定、欠定、适定方程组,是一种普适的方法;利用向量空间概念求解线性方程组,更能从本质上把握线性方程组的解的性质。应用MATL B语言编程可以轻松实现这些求解方法。实用共轭梯度法主要有以下优点 :1)算法中,系数矩阵A的作用仅仅是用来由已知向量P产生向量 w=Ap,这就充分利用了A的稀疏性。而且,对某些提供矩阵A较为困难而由已知向量p产生向量w=Ap又十分方便的应用问题是很有益的 。2)不需要预先估计任何参数就可以计算,这与SOR方法等不同。3)每次迭代所需的计算,主要是向量之间的计算,便于并行运算5 附录MATLAB语言是一种以矩阵运算为基础的计算语言。对于实现线方程组的求解非常方便。对一个四元一次方程组的求解.可以用克拉默法则和逆阵乘积法来实现,程序如下 : tic;D=[1 1 1 1;1 2 -1 4;2 -3 -1 -5;3 1 2 11]det(D)b=[5-2-2 0];Dl=[5 1 l l;-2 2 -1 4;-2 -3 -1 -5;0 1 2 11];D2=[1 5 1 1 -2 -1 4;2 -2 -1 -5;3 0 2 11];D3=[1 1 5 1;1 2-2 4;2 -3 -2 -5;3 1 0 11];D4=[1 1 1 5;1 2 -1 -2;2 -3 -1 -2;3 l 2 0];X1= det (D1)/det(D);X2= det (D2)/det(D);X3= det (D3)/det(D);X4= det (D4)/det(D);X=inv(D)*b;Toc其中克拉默法则用行列式除法( Xi=det(Di)/det(D))来实现;逆阵乘积法用(X=inv(D)*b)来实现;det(D)是系数矩阵D的行列式运算;inv(D)是D的逆阵运算。上例中.系数矩阵D不为零,可以用克拉默法则和逆阵乘积法来求解。当系数行列式为零时,只能用初等变换来求解。对于初等变换.利用阶梯生成函数命rreff也可以轻松地实现,举例如下 :tic;A=[3-4 3 2 -1;0
-6 0 -3 -3;4 -3 4 2 -2;1 1 1 0 -1;-2 6 -2 1 3];b=[2 -3 2 0 1];B=[A,b],[UB,ip]=rref(B)U0=UB([1:5],[1:5]),d=UB(:,6)tocUB为经过初等变换以后的行阶梯矩阵.可以轻松地求出方程组的解.tic、toc为计时函数,CPU-pentium1.7GHz、512MB内存、MATLAB6.5的运算环境中,两个程序的运行时间分别为0.031s和0.015s。可见,MATLAB语言实现线性方程组的求解具有程序简单、直观的特 点。同时还具有计算效率高的优点。在实际计算中摆脱了系数矩阵阶数、未知元数等的限制 。参考文献[1]同济大学应用数学系.数学——线 性代数[M].4版.北京 :高等教育出版社.2003.[2]陈怀琛,龚杰民.线性代数实践及MATLABA人民[M].北京:电子工业出版社.2005.
范文八:Hilbert矩阵Hn=[hi,j]∈Rn*n,其元素hi,j=1/(i+j-1),分别对n=2,3,4…... (1)计算cond(Hn)∞,分析条件数作为n的函数如何让变化?(2)令x=(1,1…..1)T∈Rn,计算bn=Hn*x,然后用Gauss消去法或cholesky方法解方程组bn=Hn*x,解出的解为x,
计算剩余变量rn=bn - Hn*x和误差向量;▽x=x-x。(3)分析对每个n,x分量的有效数字如何随n变化,此变化与条件数有何联系。当n为多大时绝对误差达到100%(即X连一位有效数字也没有了)。解:希尔伯特矩阵是著名的病态矩阵。n阶Hilbert矩阵为1/2?1?1/21/3?Hn??????1/n1/(n?1)其条件数 Cond(Hn)∞ 如下表所示?1/n??1/(n?1)?? ?????1/(2n?1)?Input(‘input n=’); H=hilb(n) Cond(H,inf)随着n的增大,矩阵条件数迅速增加。猜测:希尔伯特矩阵条件数以指数规律增长。即,设矩阵阶数为n,有Cond(Hn) ∞≈ exp( a n + b )用数据拟合的方法验证。问题分析:由于选择拟合函数为指数函数,直接列出超定方程组将是非线性的方程组。为了便于计算,对表中的条件数做对数变换,问题转化为线性拟合问题ln[Cond(Hn) ∞] = a n + b,( n = 2,3,4,5,6)实际操作时使用MATLAB的多项式拟合命令。线性拟合图形如下线性函数的两个系数分别为a = 3.3216,b = – 3.3474故指数拟合函数为:Cond(Hn) ∞ ≈ exp(3.3216 n –3.3474 )拟合函数的残差向量为MATLAB程序段如下 C=[]; for k=2:6H=hilb(k);C=[C,cond(H,inf)]; endLC=log(C);n=2:6; P=polyfit(n,LC,1)plot(n,LC,'o',n,polyval(P,n)) residure=C-exp(polyval(P,n))2. x=(1,1…..1)T∈Rn,.,令b=Hnx,用高斯消去法求解方程组Hnx = b解出x。。函数定义部分function x=gauss(A,b)%x=gauss(A,b)n=length(A);a=[A,b];for k=1:n-1maxa=max(abs(a(k:n,k)));if maxa==0
endfor i=k:nif abs(a(i,k))==maxay=a(i,k:n+1);a(i,k:n+1)=a(k,k:n+1);a(k,k:n+1)=y;
endfor i=k+1:nl(i,k)=a(i,k)/a(k,k);a(i,k+1:n+1)=a(i,k+1:n+1)-l(i,k).*a(k,k+1:n+1);
end%回代if a(n,n)==0returnendx(n)=a(n,n+1)/a(n,n);for i=n-1:-1:1x(i)=(a(i,n+1)-sum(a(i,i+1:n).*x(i+1:n)))/a(i,i);
endn=input('input n='); >> E=ones(n,1); >> H=hilb(n); >> b=H*E;>> x=gauss(H,b)err=norm(b-(x*H)',inf) err=norm(x-E’,inf) err=norm((x-E’)/E,inf)
范文九:线性方程组的数值解法及其应用一、问题描述现实中的问题大多数是连续的,例如工程中求解结构受力后的变形,空气动力学中计算机翼周围的流场,气象预报中计算大气的流动。这些现象大多是用若干个微分方程描述。用数值方法求解微分方程(组),不论是差分方法还是有限元方法, 通常都是通过对微分方程(连续的问题, 未知数的维数是无限的)进行离散,得到线性方程组(离散问题, 因为未知数的维数是有限的)。因此线性方程组的求解在科学与工程中的应用非常广泛。经典的求解线性方程组的方法一般分为两类:直接法和迭代法。二、基本要求1)掌握用MATLAB软件求线性方程初值问题数值解的方法;2)通过实例学习用线性方程组模型解决简化的实际问题;3)了解用高斯赛德尔列主元消去法和雅可比迭代法解线性方程组。三、测试数据1) 直接法:A=[0.002 52.88;4.573 -7.290];b=[52.90;38.44];2) 迭代法:A=[10 -1 -2;-1 10 -2;-1 -1 5];b=[7.2;8.3;4.2];四、算法程序及结果1)function[RA,RB,n,x]=liezy1(A,b)B=[A b];n=length(b);RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0,disp('因为RA~=RB,所以此方程组无解.')
return end if RA==RB
if 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);
endb=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);
elsedisp('因为RA=RB测试:A=[0.002 52.88;4.573 -7.290]; >> b=[52.90;38.44]; >> [RA,RB,n,x]=liezy1(A,b) 因为RA=RB=n,所以此方程组有唯一解. RA =2RB =
2x =10.00001.00002)function Jacobi(A,b,x0,P,error,max1) [n n]=size(A);x=zeros(n,1); for k=1:max1for j=1;nx(j)=(b(j)-A(j,[1:j-1,j+1:n])*x0([1:j-1,j+1:n]))/A(j,j);
xerrx=norm(x-x0,P);
x0=x;x1=A\b;
if(errxdisp('迭代次数k,精确解x1和近似解x分别是:')
end endif(errx>=error)disp('请注意:Jacobi迭代次数已经超过最大迭代次数max1.') end
测试:A=[10 -1 -2;-1 10 -2;-1 -1 5]; >>b=[7.2;8.3;4.2]; >>x0=[0;0;0];>>Jacobi(A,b,x0,inf,0.001,100)
3 x =0.7200
0迭代次数k,精确解x1和近似解x分别是:
k =2x1 =1.1000
0五、应用举例1)营养学家配制一种具有1200卡,30g蛋白质及300mg维生素C的配餐。有3种食物可供计算所需果冻、鲜鱼、牛肉的数量。 设果冻、鲜鱼、牛肉的数量为x1,x2,x3则 20x1+100x2+200x3=x2+2x3=30 30x1+20x2+10x3=3002)交通流量问题下图给出了某城市部分单行街道的交通流量(每小时通过的车辆数)图中有6个路口,已有9条街道记录了当天的平均车流量。另有7处的平均车流量未知,试利用每个路口的进出车流量相等关系推算这7处的平均车流量。其数学模型为线性方程组:3)闭合经济问题一个木工,一个电工,一个油漆工,三人相互同意彼此装修他们自己的房子。在装修之前,他们达成了如下协议:(1)每人总共工作十天(包括给自已家干活在内); (2)每人的日工资根据一般的市价在60元~80元之间;(3)每人的日工资数应使得每人的总收入与总支出相等。下面的表是他们协商后制定出的工作天数的分配方案:设木工的日工资为x1,电工的日工资为x2,油漆工的日工资为x3。则木工的10个工作日总收入应该为 10x1 ,而木工、电工及油漆工三人在木工家工作的天数分别为:2天, 1天, 6天,按日工资累计木工的总支出为 2x1 + x2 + 6x3。于是木工的收支平衡可描述为等式同理,可建立描述电工,油漆工各自的收支平衡关系的另外两个等式,将三个等式联立,可得描述实际问题的方程组。六、实训小结通过本次实训,我们明白了如何用MATLAB解线性方程组,也了解了线性方程组的两种解法:直接法和迭代法,如高斯列主元消去法、雅克比迭代法,在实训过程中也碰到了一些困难,但我们都一一解决了,既增加了知识面,也加强了理解,此次实训我们分工明确,收获颇丰,更加明白了线性方程组,或者说数学在生活中的重要应用。七、参考文献数值分析 中国科技大学出版 朱晓琳编数学软件 MATLAB 5.3
范文十:实验六:线性方程组迭代解法1)实验目的o 熟悉Matlab编程;o 学习线性方程组迭代解法的程序设计算法 2)实验题目1.研究解线性方程组Ax=b迭代法收敛速度。A为20阶五对角距阵??3?1???2??1?4A???????????12312?1?4?141?2?3?1?21?4141?2??3121?4?14?1?2?3?12????????? 1???4?1???2?3??要求:(1)选取不同的初始向量x0 及右端向量b,给定迭代误差要求,用雅可比迭代和高斯-赛德尔迭代法求解,观察得到的序列是否收敛?若收敛,记录迭代次数,分析计算结果并得出你的结论。(2)用SOR迭代法求解上述方程组,松弛系数ω取12.给出线性方程组Hnx?b,其中系数矩阵Hn为希尔伯特矩阵:Hn??hij???n?n,hij?Ti,i,j?1,2,?,n.i?j?1假设x*??1,1,?,1???n,b?Hnx.若取n?6,8,10,分别用雅可比迭代法及SOR迭代(??1,1.25,1.5)求解,比较计算结果。 3)实验原理与理论基础1.雅克比(Jacobi)迭代法算法设计:①输入矩阵a与右端向量b及初值x(1,i);
②按公式计算得xi(k?1)1?aii??n?(k)?b?ax?i?ijj?j?1??j?i??(i?1,2,?,n)2.高斯――赛得尔迭代法算法设计:1. 输入矩阵a与右端向量b及初值x(1,i).2.x(k?1)i1?aiii?1??bi??aijx(jk?1)??j?1?(k)??ax?ijj? (i = 1, 2,…, n)j?i?1?n3.超松驰法算法设计:①输入矩阵a与右端向量b及初值x(1,i)。②x(k?1)i?(1??)x(k)i?bi??aijxaii?j?1????i?1(k?1)j?j?i?1?anijx(k)j??,0???2 ??4)实验内容第一题实验程序: 1.雅克比迭代法: function []=yakebi(e) %输入矩阵a与右端向量b。 for i=1:20
a(i,i)=3; endfor i=3:20
for j=i-2a(i,j)=-1/4;
a(j,i)=-1/4;
end endfor i=2:20
for j=i-1a(i,j)=-1/2;
a(j,i)=-1/2;
end endb=[2.2 1.7 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.7 2.2]; k=1;n=length(a); for i=1:nx(1,i)=1;%数组中没有第0行。 endwhile k>=1
m=0;%此步也可以用ifj~=i条件判定一下。
for j=1:(i-1)m=m+a(i,j)*x(k,j);
endfor j=(i+1):nm=m+a(i,j)*x(k,j);
endx(k+1,i)=(b(i)-m)/a(i,i);
l=0;%判定满足条件使循环停止迭代。
for i=1:nl=l+abs(x(k+1,i)-x(k,i));
if l%输出所有的x的值。
x(k+1,:)k2.高斯—赛德尔迭代法: function []=gaoshisaideer(e) for i=1:20
a(i,i)=3; endfor i=3:20
for j=i-2a(i,j)=-1/4;
a(j,i)=-1/4;
end endfor i=2:20
for j=i-1a(i,j)=-1/2;
a(j,i)=-1/2;
end endb=[2.2 1.7 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.7 2.2]; k=1;n=length(a); for i=1:nx(1,i)=0;%数组中没有第0行。 endwhile k>=1
p=0;q=0;for j=1:(i-1)p=p+a(i,j)*x(k+1,j);
endfor j=(i+1):nq=q+a(i,j)*x(k,j);
endx(k+1,i)=(b(i)-q-p)/a(i,i);
l=0;%判定满足条件使循环停止迭代。
for i=1:nl=l+abs(x(k+1,i)-x(k,i));
if l%输出所有的x的值。
x(k+1,:) k3.SOR迭代法程序:function []=caosongci(e,w) for i=1:20
a(i,i)=3; endfor i=3:20
for j=i-2a(i,j)=-1/4;
a(j,i)=-1/4;
end endfor i=2:20
for j=i-1a(i,j)=-1/2;
a(j,i)=-1/2;
end endb=[2.2 1.7 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.7 2.2]; k=1;n=length(a); for i=1:nx(1,i)=0;%数组中没有第0行。 endwhile k>=1if w>=2||w'请重新输入w的值,w在1与2之间';break
p=0;q=0;for j=1:(i-1)p=p+a(i,j)*x(k+1,j);
endfor j=i:nq=q+a(i,j)*x(k,j);
endx(k+1,i)=x(k,i)+w*(b(i)-q-p)/a(i,i);
l=0;%判定满足条件使循环停止迭代。
for i=1:nl=l+abs(x(k+1,i)-x(k,i));
if l%输出所有的x的值。
x(k+1,:) k第二题实验程序: 1.雅克比迭代法:function X = p211_1_JJ(n)Hn = GET_Hn(n); b = GET_b(n); temp = 0;X0 = zeros(1, n); X_old = zeros(1, n); X_new = zeros(1, n);disp('Now Jacobi method!');disp('Start with the vector that (0, 0, 0, ...)^T'); for i = 1:nfor k = 1:nX_old = X_
for j = 1:n
if(j ~= i)temp = temp + Hn(i, j) * X_old(j);
endendX_new(i) = (b(i) - temp) / Hn(i, i);
end endX = X_ end2.SOR迭代法:function X = p211_1_SOR(n, w) Hn = GET_Hn(n); b = GET_b(n); temp01 = 0; temp02 = 0;X0 = zeros(1, n); X_old = zeros(1, n); X_new = zeros(1, n);disp('Now Successive Over Relaxtion method!');disp('Start with the vector that (0, 0, 0, ...)^T'); for i = 1:n
for k = 1:nX_old = X_
temp01 = 0;
temp02 = 0;
for j = 1:n
if(jtemp01 = temp01 + Hn(i, j) * X_new(j);
endif(j > i)temp02 = temp02 + Hn(i, j) * X_old(j);
endX_new(i) = w * (b(i) - temp01 -
temp02) / Hn(i, i) + X_old(i); endX = X_ end5)实验结果 第一题实验结果: 1.雅克比迭代法: 输入:>> b=[2.2 1.7 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.7 2.2]; yakebi(0.00001) 结果: ans =Columns 1 through 120.9793
1.0000Columns 13 through 200.9999
122.高斯—赛德尔迭代法: 此时初值全取1; 输入:>> b=[2.2 1.7 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.7 2.2]; >> gaoshisaideer(0.00001) 结果:ans =Columns 1 through 120.9793
Columns 13 through 200.9999
0.9793 k =14此时初值全取1;输入:>> b=[2.5 1.9 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.9 2.5]; gaoshisaideer(0.00001) 结果: ans =Columns 1 through 121.0969
1.0001Columns 13 through 201.0003
143.SOR迭代法:>> caosongci(0.) ans =Columns 1 through 120.9793
1.0000Columns 13 through 200.9999
0.97870.9793k =
11>> caosongci(0.) ans =Columns 1 through 120.9793
1.0000Columns 13 through 200.9999
12>> caosongci(0.) ans =Columns 1 through 120.9793
Columns 13 through 200.9999
15>> caosongci(0.) ans =Columns 1 through 120.9793
Columns 13 through 200.9999
19>> caosongci(0.) ans =Columns 1 through 120.9793
Columns 13 through 200.9999
25>> caosongci(0.)0.9989
1.00000.9989
1.00000.9989
1.00000.9989
0.9787ans =Columns 1 through 120.9793
1.0000Columns 13 through 200.9999
34>> caosongci(0.) ans =Columns 1 through 120.9793
1.0000Columns 13 through 200.9999
47>> caosongci(0.) ans =Columns 1 through 120.9793
1.0000Columns 13 through 200.9999
73>> caosongci(0.) ans =Columns 1 through 120.9793
1.0000Columns 13 through 200.9999
150第二题实验结果: 1.雅克比迭代法:
>> p211_1_JJ(6)Now Jacobi method!Start with the vector that (0, 0, 0, ...)^T0.9995
0.9787ans =2.4500
>> p211_1_JJ(8)
Now Jacobi method!Start with the vector that (0, 0, 0, ...)^T ans =2.7179
0.1995>> p211_1_JJ(10)Now Jacobi method!Start with the vector that (0, 0, 0, ...)^T ans =Columns 1 through 92.9290
Column 100.19512.SOR迭代法:n=6, ω=1,1.25,1.5的时候 >> p211_1_SOR(6, 1)Now Successive Over Relaxtion method!Start with the vector that (0, 0, 0, ...)^T ans =2.4500
0.2071 >> p211_1_SOR(6, 1.25)Now Successive Over Relaxtion method!Start with the vector that (0, 0, 0, ...)^T ans =3.0625
0.2097 >> p211_1_SOR(6, 1.5)Now Successive Over Relaxtion method!Start with the vector that (0, 0, 0, ...)^T ans =3.6750
-0.0384 与n=8, ω=1,1.25,1.5的时候 >> p211_1_SOR(8, 1)Now Successive Over Relaxtion method!Start with the vector that (0, 0, 0, ...)^T ans =2.7179
0.1995>> p211_1_SOR(8, 1.25)Now Successive Over Relaxtion method!Start with the vector that (0, 0, 0, ...)^Tans =3.3973
0.2042>> p211_1_SOR(8, 1.5)Now Successive Over Relaxtion method!Start with the vector that (0, 0, 0, ...)^Tans =4.0768
0.1275与n=10, ω=1,1.25,1.5的时候>> p211_1_SOR(10, 1)Now Successive Over Relaxtion method!Start with the vector that (0, 0, 0, ...)^Tans =Columns 1 through 92.9290
0.2325Column 100.1951>> p211_1_SOR(10, 1.25)Now Successive Over Relaxtion method!Start with the vector that (0, 0, 0, ...)^Tans =Columns 1 through 93.6612
0.2363Column 100.1984>> p211_1_SOR(10, 1.5)Now Successive Over Relaxtion method!Start with the vector that (0, 0, 0, ...)^Tans =Columns 1 through 94.3935
0.2819Column 100.17666)实验结果分析与小结本次实习主要是学会应用雅克比迭代法、高斯—赛德尔迭代法、SOR迭代法三种迭代法,并且了解三种迭代法的性质以及迭代精度等。第一题中取的b对于雅克比迭代法、高斯――赛得尔迭代法都是收敛的,对于相同的初值与右端向量明显可以看出高斯――赛得尔迭代法比雅克比迭代法快。由第二题可得出对于SOR迭代方法选择不同的松弛因子,收敛次数大大不同,而当松弛因子为1.1时,在同等条件下迭代最快。}

我要回帖

更多关于 pcb覆铜厚度 的文章

更多推荐

版权声明:文章内容来源于网络,版权归原作者所有,如有侵权请点击这里与我们联系,我们将及时删除。

点击添加站长微信