解方程组基本思想.docx
- 文档编号:8156365
- 上传时间:2023-01-29
- 格式:DOCX
- 页数:32
- 大小:105.68KB
解方程组基本思想.docx
《解方程组基本思想.docx》由会员分享,可在线阅读,更多相关《解方程组基本思想.docx(32页珍藏版)》请在冰豆网上搜索。
解方程组基本思想
四:
基本方法
基本思路将在解题的过程中得到体现。
1.(求线性方程组的唯一解或特解),这类问题的求法分为两类:
一类主要用于解低阶稠
密矩阵——直接法;一类是解大型稀疏矩阵——迭代法。
1.1利用矩阵除法求线性方程组的特解(或一个解)
方程:
AX=b,解法:
X=A\b,(注意此处’\’不是’/’)
例1-1 求方程组的解。
解:
A=;=;b=(1,0,0,0,1)’
由于>>rank(A)=5,rank()=5 %求秩,此为R(A)=R()>=n的情形,有唯一解。
>>X=A\b %求解X=(2.2662,-1.7218,1.0571,-0.5940,0.3188)’或用函数rref
求解,>>sv=rref(A:
b);所得sv的最后一列即为所要求的解。
1.2利用矩阵的LU、QR和cholesky分解求方程组的解
这三种分解,在求解大型方程组时很有用。
其优点是运算速度快、可以节省磁盘空间、节省内存。
I)LU分解又称Gauss消去分解,可把任意方阵分解为下三角矩阵的基本变换形式(行交换
)和上三角矩阵的乘积。
即A=LU,L为下三角阵,U为上三角阵。
则:
A*X=b 变成L*U*X=b
所以X=U\(L\b) 这样可以大大提高运算速度。
命令 [L,U]=lu(A)
在matlab中可以编如下通用m文件:
在Matlab中建立M文件如下
%exp1.m
A;b;
[L,U]=lu(A);
X=U\(L\b)
II)Cholesky分解
若A为对称正定矩阵,则Cholesky分解可将矩阵A分解成上三角矩阵和其转置的乘积,即:
其中R为上三角阵。
方程 A*X=b 变成 所以
在Matlab中建立M文件如下
%exp2.m
A;b;
[R’,R]=chol(A);
X=R\(R’\b)
III)QR分解
对于任何长方矩阵A,都可以进行QR分解,其中Q为正交矩阵,R为上三角矩阵的初等变换形
式,即:
A=QR
方程 A*X=b 变形成 QRX=b
所以 X=R\(Q\b)
上例中 [Q,R]=qr(A)
X=R\(Q\B)
在Matlab中建立M文件如下
%exp3.m
A;b;
[Q,R]=qr(A);
X=R\(Q\b)
2.求线性齐次方程组的通解(A*X=0)
在Matlab中,函数null用来求解零空间,即满足A•X=0的解空间,实际上是求出解空
间的一组基(基础解系)。
在Matlab中建立M文件如下
%exp4.m
formatrat %指定有理式格式输出
A;b=0;
r=rank(A);
bs=null(A,‘r’);%一组基含(n-r)个列向量
%k,k,……,k
%X=k*bs(:
1)+k*bs(:
2)+……+k*bs(:
n-r)方程组的通解
pretty(X) %让通解表达式更加精美
3 求非齐次线性方程组的通解(A*X=b)
非齐次线性方程组需要先判断方程组是否有解,若有解,再去求通解。
因此,步骤为:
第一步:
判断AX=b是否有解,(利用基本思路的第一条)
若有解则进行第二步
第二步:
求AX=b的一个特解
第三步:
求AX=0的通解
第四步:
AX=b的通解为:
AX=0的通解加上AX=b的一个特解。
在Matlab中建立M文件如下
%exp4.m
clearall
A;b; %输入矩阵A,b
[m,n]=size(A);
R=rank(A);
B=[Ab];
Rr=rank(B);
formatrat
ifR==Rr&R==n %n为未知数的个数,判断是否有唯一解
x=A\b;
elseifR==Rr&R x=A\b %求特解 C=null(A,r) %求AX=0的基础解系,所得C为n-R列矩阵,这n-R列即为对 %应的基础解系 %这种情形方程组通解xx=k(p)*C(: P)(p=1…n-R) elseX=Nosolution! %判断是否无解 end 第3章线性方程组的迭代解法 3.1实验目的 理解线性方程组计算机解法中的迭代解法的求解过程和特点,学习科学计算的方法和简单的编程技术。 3.2 概念与结论 1. n阶线性方程组 如果未知量的个数为n,而且关于这些未知量x1,x2,…,xn的幂次都是一次的(线性的)那末,n个方程 a11x1+a12x2+…+a1nxn=b1 ┆┆┆ (1) an1x1+an2x2+…+annxn=bn 构成一个含n个未知量的线性方程组,称为n阶线性方程组。 其中,系数a11,…,a1n,a21,…,a2n,…,an1,…,ann和b1,…,bn都是给定的常数。 方程组 (1)也常用矩阵的形式表示,写为 Ax=b 其中,A是由系数按次序排列构成的一个n阶矩阵,称为方程组的系数矩阵,x和b都是n维向量,b称为方程组的右端向量。 2.n阶线性方程组的解 使方程组 (1)中每一个方程都成立的一组数x1*,x2*,…,xn*称为式 (1)的解,把它记为向量的形式,称为解向量. 3.向量范数的三种常用范数 4.矩阵的四种常用范数 5.谱半径 设nn阶矩阵A的特征值为i(i=1,2,3……n),则称 (A)=MAX|i| 1in 为矩阵A的谱半径. 矩阵范数与谱半径之间的关系为: (A)||A|| 6.严格(行)对角占优阵A 如果矩阵A=(aij)满足 n |aii|>|aij|i=1,2,……n, j=1,ji 则称方阵A是严格(行)对角占优的. 7.收敛定理 对任意初始向量x(0)及任意右端向量g,由迭代x(k+1)=Bx(k)+g产生的迭代向量序列{x(k)}收敛的充要条件是谱半径(B)<1 8.收敛判别条件 判别条件1: 若||B||<1,则迭代x(k+1)=Bx(k)+g对任何初始向量x(0)都收敛. 判别条件2: 如果A为严格对角占优阵,则其Jacobi迭代和Seidel迭代对任何初始向量x(0)都收敛。 判别条件3: 如果A为对称正定阵,则其Seidel迭代对任何初始向量x(0)都收敛。 9.迭代法的误差估计 若||B||<1,则对迭代格式x(k+1)=Bx(k)+g 有 3.3程序中Mathematica语句解释 a*matrix数a与矩阵matrix相乘 matrix1+matrix2矩阵matrix1和矩阵matrix2相加(注意矩阵的大小相同) matrix1.matrix2矩阵matrix1和矩阵matrix2相乘(注意矩阵乘法的规则) Transpose[matrix]求矩阵matrix转置 Inverse[matrix]求矩阵(方阵)matrix的逆 DiagonalMatrix[list]使用列表list中的元素生成一个对角矩阵. IdentityMatrix[n]生成n阶单位矩阵 Max[x]求向量x中元素的最大值 3.4方法、程序、实验 解线性方程组的迭代法是将线性方程组Ax=b化为等价线性方程组 x=Bx+f 再由矩阵迭代格式 x(k+1)=Bx(k)+f 构造向量序列{x(k)}来求线性方程组解的。 如果得出的向量序列{x(k)}收敛至某个向量x*,则可得该向量x*就是所求方程组Ax=b的准确解. 线性方程组的迭代法主要有Jocobi迭代法、Seidel迭代法和超松弛(Sor)迭代法。 1.Jocobi迭代法 1)Jocobi迭代法的构造过程 假设aii0,依次在第i个方程解出xi,i=1,2,,n并令 cij=-aij/aii(ij),gi=bi/aii 就得到如下Jocobi迭代格式: x1(k+1)=c12x2(k)+c13x3(k)++c1nxn(k)+g1 x2(k+1)=c21x1(k)+c23x3(k)++c2nxn(k)+g2 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 xn(k+1)=cn1x1(k)+cn2x2(k)++cn(n-1)xn-1(k)+gn 若令 则有Jocobi迭代的矩阵格式: x(k+1)=BJx(k)+gJ BJ称为Jocobi迭代矩阵。 Jocobi迭代可以写成如下紧凑格式: 在给定初始迭代向量x(0)后就可以进行Jocobi迭代求解了。 2)Jacobi迭代算法 1.输入变量个数n、初值向量x(0)、迭代精度eps、系数矩阵A、常数项b和迭代最大次数nmax 2Fori=1,2,…,n 2.1如果|aii| 3.BjE-D-1A 4.gjD-1b 5.Fork=1,2,…,nmax 5.1xBj.x0+gj 5.2如果||x-x0|| 6.如果||x-x0||>eps,输出迭代失败,终止。 3)Jacobi迭代法程序 Clear[a,b,x]; nmax=500; n=Input[“线性方程组阶数n=”]; a=Input["系数矩阵A="]; b=Input["常数项b="]; x0=Input["输入迭代初值向量x0"]; eps1=0.000001; eps=Input["输入精度控制eps="]; Do[If[Abs[a[[i,i]]] If[t1==1, Print["Jacobi迭代法失效"], d=DiagonalMatrix[Table[a[[i,i]],{i,1,n}]]; d1=Inverse[d]; bj=IdentityMatrix[n]-d1.a; gj=d1.b; Do[x=bj.x0+gj; err=Max[Abs[x-x0]]; Print["x=",x//N,"i=",i,"err=",err//N]; If[N[err] {i,1,nmax}]; If[err>=eps,Print["迭代失败"]] ] 说明本程序用于求线性方程组Ax=b的解。 程序执行后,先通过键盘输入线性方程组阶数n、系数矩阵A、常数项b、迭代初值向量x0和输入精度控制eps,程序即可给出每次迭代的次数和对应的迭代向量序列x(k),其中最后输出的结果即为所求的根。 如果迭代超出500次还没有求出满足精度的根则输出迭代失败提示,如果出现主对角线元素aii=0给出Jacobi迭代法失效提示。 程序中变量说明 x0: 存放初始向量和迭代过程中的向量x(k) x: 存放迭代过程中的向量x(k+1) nmax: 存放迭代允许的最大次数 err: 存放误差||x-x0|| t1: 临时变量 注: 迭代最大次数可以修改为其他数字。 4)例题与实验 例1.用Jacobi迭代法解如下线性方程组 5x1+2x2+x3=-12 -x1+4x2+2x3=20 2x1-3x2+10x3=3 要求误差||x(k+1)-x(k)||<10-4,并用取不同初值的方法实验观察迭代收敛的情况。 解: 执行Jacobi迭代法程序后在输入的四个窗口中按提示分别输入: 3、{{5,2,1},{-1,4,2},{2,-3,10}}、{-12,20,3}、{0,0,0}、0.0001 每次输入后用鼠标点击窗口的“OK”按扭,得如下输出结果: x={-2.4,5.,0.3}i=1err=5. x={-4.46,4.25,2.28}i=2err=2.06 x={-4.556,2.745,2.467}i=3err=1.505 x={-3.9914,2.6275,2.0347}i=4err=0.5646 x={-3.85794,2.9848,1.88653}i=5err=0.3573 x={-3.97123,3.09225,1.96703}i=6err=0.113286 x={-4.03031,3.02368,2.02192}i=7err=0.0685705 x={-4.01386,2.98146,2.01316}i=8err=0.042216 x={-3.99522,2.98995,1.99721}i=9err=0.0186374 x={-3.99542,3.00259,1.99603}i=10err=0.0126367 x={-4.00024,3.00313,1.99986}i=11err=0.0048186 x={-4.00122,3.00001,2.00099}i=12err=0.00312067 x={-4.0002,2.9992,2.00025}i=13err=0.00102319 x={-3.99973,2.99983,1.9998}i=14err=0.000625697 x={-3.99989,3.00017,1.99989}i=15err=0.00034136 x={-4.00005,3.00008,2.00003}i=16err=0.000155236 x={-4.00004,2.99997,2.00003}i=17err=0.000106099 x={-4.,2.99997,2.}i=18err=0.0000414468 此结果说明迭代18次,求得误差为err=0.0000414468的近似解,最后显示的近似解向量为x={-4.,2.99997,2.},它表示所求解为 x1=-4,x2=2.99997,x3=2。 本题的准确解为 x1=-4,x2=3,x3=2。 如果将如上输入的初值改为{21,-18,30},执行Jacobi迭代法程序后得如下输出结果: x={-1.2,-4.75,-9.3}i=1err=39.3 x={1.36,9.35,-0.885}i=2err=14.1 x={-5.963,5.7825,2.833}i=3err=7.323 x={-5.2796,2.09275,3.22735}i=4err=3.68975 x={-3.88257,2.06642,1.98375}i=5err=1.39703 x={-3.62332,3.03749,1.69644}i=6err=0.97106 x={-3.95428,3.24595,1.93591}i=7err=0.330963 x={-4.08556,3.04347,2.06464}i=8err=0.202475 x={-4.03032,2.94629,2.03015}i=9err=0.0971858 x={-3.98455,2.97734,1.98995}i=10err=0.0457716 x={-3.98893,3.00889,1.99011}i=11err=0.0315451 x={-4.00158,3.00771,2.00045}i=12err=0.0126504 x={-4.00318,2.99938,2.00263}i=13err=0.00833246 x={-4.00028,2.99789,2.00045}i=14err=0.00289753 x={-3.99925,2.99971,1.99942}i=15err=0.0018145 x={-3.99977,3.00048,1.99976}i=16err=0.000770763 x={-4.00014,3.00018,2.0001}i=17err=0.000375926 x={-4.00009,2.99992,2.00008}i=18err=0.000261658 x={-3.99998,2.99994,1.99999}i=19err=0.000107579 x={-3.99997,3.00001,1.99998}i=20err=0.0000714045 从计算结果可以看到虽然所取的初值不同,但迭代总是收敛相同的结果。 考察本题系数矩阵可以发现它是严格对角占优矩阵,由收敛判别法,Jacobi迭代对任意初值都是收敛的,我们通过实际计算验证了这个结果。 例2.通过Jacobi迭代序列观察用Jacobi迭代法解如下线性方程组 x1+2x2+2x3=-12 -x1+4x2+x3=20 2x1-3x2+x3=3 的收敛性。 解: 执行Jacobi迭代法程序后在输入的四个窗口中按提示分别输入: 3、{{1,2,2},{-1,4,1},{2,-3,1}}、{-12,20,3}、{0,0,0}、0.0001 每次输入后用鼠标点击窗口的“OK”按扭,得如下输出结果: x={-12.,5.,3.}i=1err=12. x={-28.,1.25,42.}i=2err=39. x={-98.5,-12.5,62.75}i=3err=70.5 x={-112.5,-35.3125,162.5}i=4err=99.75 ……………………………………………. {246.719,-120.176,696.719}i=8err=763.5 x={-1165.09,-107.5,-850.965}i=9err=1547.68 x={1904.93,-73.5303,2010.67}i=10err=3070.02 x={-3886.28,-21.4355,-4027.45}i=11err=6038.12 x={8085.77,40.2917,7711.26}i=12err=11972.1 x={-15515.1,98.6279,-16047.7}i=13err=23758.9 x={31886.1,138.141,31329.1}i=14err=47401.2 x={-62946.5,144.247,-63354.7}i=15err=94832.5 x={126409.,107.068,126329.}i=16err=189683. x={-252883.,25.0775,-252494.}i=17err=379292. x={504925.,-92.4305,505845.}i=18err=758339. ……………………………………. 通过观察迭代过程中的误差是不断变大的特点可以知道本题的Jacobi迭代序列是不收敛的,因此,本题线性方程组不能用Jacobi迭代法求解。 这个实验说明并不是每个线性方程组都能用Jacobi迭代法求解。 2. Seidel迭代 1)Seidel迭代的构造过程 为了加快收敛速度,同时节省计算机的内存,对Jocobi迭代作如下的改进: 每算出一个分量的近似值,立即用到下一个分量的计算中去,即用迭代格式: x1(k+1)=c12x2(k)+c13x3(k)++c1nxn(k)+g1 x2(k+1)=c21x1(k+1)+c23x3(k)++c2nxn(k)+g2 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 。 xn(k+1)=cn1x1(k+1)+cn2x2(k+1)++cn(n-1)xn-1(k+1)+gn 如上迭代可以写成如下紧凑格式: 这样所得的迭代法就称为Seidel迭代法。 利用Ax=b及A=L+D+U,其中D为对角矩阵,L,U分别为严格下,上三角矩阵.则有Seidel迭代法的矩阵形式为 x(k+1)=Bsx(k)+gs其中: Bs=-(L+D)-1U,gs=D-1b Seidel迭代矩阵为Bs=-(L+D)-1U。 在给定初始迭代向量x(0)后就可以进行Seidel迭代求解了。 Jacobi迭代和Seidel迭代格式可表述为统一形式: 2)Seidel迭代算法 1.输入变量个数n、初值向量x(0)、迭代精度eps、系数矩阵A、常数项b和迭代最大次数nmax 2.Fori=1,2,…,n 2.1如果|aii| 3.Fori=1,2,…,nmax 3.1Fori=1,2,…,n 3.2如果||x(k+1)-x(k)|| 4.如果||x-x0|| 3)Seidel迭代法程序 Clear[a,b,x]; nmax=500; n=Input[“线性方程组阶数n=”]; a=Input["系数矩阵A="]; b=Input["常数项b="]; x0=Input["输入迭代初值向量x0"]; eps1=0.000001; eps=Input["输入精度控制eps="]; x=x0; Do[If[Abs[a[[i,i]]] If[t1==1, Print["Seidel迭代法失效"], Do[Do[u1=Sum[a[[i,j]]*x[[j]],{j,1,i-1}]; u2=Sum[a[[i,j]]*x0[[j]],{j,i+1,n}]; x[[i]]=(b[[i]]-u1-u2)/a[[i,i]], {i,1,n}]; err=Max[Abs[x-x0]]; Print["x=",x//N,"k=",k,"err=",err//N]; If[N[err] {k,1,50}]; If[err>=eps,Print["迭代失败"]] ] 说明本程序用于求线性方程组Ax=b的解。 程序执行后,先通过键盘输入线性方程组阶数n、系数矩阵A、常数项b、迭代初值向量x0和输入精度控制eps,程序即可给出每次迭代的次数和对应的迭代向量序列x(K),其中最后输出的结果即为所求的根。 如果迭代超出500次还没有求出满足精度的根则输出迭代失败提示,如果出现主对角线元素aii=0给出Jacobi迭代法失效提示。 程序中变量说明 x0: 存放初始向量和迭代
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 方程组 基本 思想
![提示](https://static.bdocx.com/images/bang_tan.gif)