数值分析上机试验——解线性方程组讲解VIP免费

实验报告课程名称数值分析实验项目名称解线性方程组实验类型上机实验学时4班级20111131学号2011113130姓名张振指导教师沈艳实验室名称理学楼407实验时间2013.12.9实验成绩预习部分实验过程表现实验报告部分总成绩教师签字日期哈尔滨工程大学教务处制实验四解线性方程组一．解线性方程组的基本思想1．直接三角分解法：将系数矩阵A转变成等价两个矩阵L和U的乘积，其中L和U分别是下三角和上三角矩阵。当A的所有顺序主子式都不为0时，矩阵A可以分解为A=LU，且分解唯一。其中L是单位下三角矩阵，U是上三角矩阵。2．平方根法：如果矩阵A为n阶对称正定矩阵，则存在一个对角元素为正数的下三角实矩阵L，使得：A=LL^T。当限定L的对角元素为正时，这种分解是唯一的，称为平方根法（Cholesky）分解。3．追赶法：设系数矩阵为三对角矩阵1122233111000000000000000000nnnnnbcabcabAabcab则方程组Ax=f称为三对角方程组。设矩阵A非奇异，A有Crout分解A=LU，其中L为下三角矩阵，U为单位上三角矩阵，记11222331100001000000010000000100,000000000000001nnnnbLU可先依次求出L，U中的元素后，令Ux=y，先求解下三角方程组Ly=f得出y，再求解上三角方程组Ux=y。4．雅克比迭代法：首先将方程组中的系数矩阵A分解成三部分，即：A=L+D+U，如图1所示，其中D为对角阵，L为下三角矩阵，U为上三角矩阵。之后确定迭代格式，X)1(k=BX)(k+f，如图2所示，其中B称为迭代矩阵，雅克比迭代法中一般记为J。（k=0,1，......）再选取初始迭代向量X)0(，开始逐次迭代。5．超松弛迭代法（SOR）它是在GS法基础上为提高收敛速度，采用加权平均而得到的新算法。选取分裂矩阵M为带参数的下三角矩阵M＝1（D－L），其中＞0为可选择的松弛因子，一般当1<<2时称为超松弛。二.实验题目及实验目的1．(第五章习题8)用直接三角分解（杜利特尔（Doolittle）分解）求线性方程组141x+251x+361x=9，131x+241x+351x=8，121x+2x+32x=8的解。2．(第五章习题9)用追赶法解三对角方程组Ax=b，其中A=2100012100012100012100012，b=00001.3．（第五章习题10）用改进的平方根法解线性方程组131321112321xxx=6544．（第六章习题7）用SOR方法解线性方程组（分别取松弛因子ω=1.03，ω=1，ω=1.1）41x-2x=1，-1x+42x-3x=4，-2x+43x=-3.精确解x*=（21，1，-21）T.要求当)(*kxx<5×106时迭代终止，并且对每一个ω值确定迭代次数.5.（第六章习题8）用SOR方法解线性方程组（取ω=0.9）51x-22x+3x=-12，-1x+42x-23x=20，21x-32x+103x=3.要求当)()1(kkxx<104时迭代终止.6．（第六章习题9）设有线性方程组Ax=b，其中A为对称正定阵，迭代公式)()1(kkxx+ω(b-A)(kx)，k=0，1，2⋯，试证明当0<ω<2时上述迭代法收敛（其中0<(A)）.7．（第六章计算实习题1）给出线性方程组Hnx=b，其中系数矩阵Hn为希尔伯特矩阵：Hnx=(hij)Rnn,hij=11ji，i，j=1，2，⋯，n.假设x*=（1，1，⋯，1）TRn，b=Hnx*.若取n=6，8，10，分别雅克比迭代法及SOR迭代（ω=1，1.25，1.5）求解.比较计算结果.三．实验手段：指操作环境和平台:win7系统下MATLABR2009a程序语言：一种类似C语言的程序语言，但比C语言要宽松得多，非常方便。四.程序1.①直接三角分解（文件ZJsanjiao.m）functionx=ZJsanjiao(A,b)[m,n]=size(A);[lu]=lu(A);s=inv(l)*[A,b];x=ones(m,1);fori=m:-1:1h=s(i,m+1);forj=m:-1:1;ifj~=ih=h-x(j)*s(i,j);endendx(i)=h/s(i,i);end②控制台输入代码：>>A=[1/4,1/5,1/6;1/3,1/4,1/5;1/2,1,2];>>b=[9;8;8];>>x=ZJsanjiao(A,b)2.①追赶法（文件ZG_SDJ.m）functionx=ZG_SDJ(a,b,c,f)%aê????????a??%bê???????é?·?μ??a??￡???êy±èaéùò???%cê?????????·?μ??a??￡???êy±èaéùò???%fê?3￡êy??bN=length(a);b=[b,0];c=[0,c];a1=zeros(N,1);b1=zeros(N,1);y=zeros(N,1);x=zeros(N,1);a1(1)=a(1);b1(1)=b(1)/a1(1);y(1)=f(1)/a1(1);forj1=2:Na1(j1)=a(j1)-c(j1)*b1(j1-1);b1(j1)=b(j1)/a1(j1);temp1=f(j1)-c(j1)*y(j1-1);y(j1)=temp1/a1(j1);endj1=N;x(j1)=y(j1);forj1=N-1:-1:1x(j1)=y(j1)-b1(j1)*x(j1+1);end②控制台输入代码：>>a=[22222];>>b=[-1-1-1-1];>>c=[-1-1-1-...