预览加载中,请您耐心等待几秒...
1/10
2/10
3/10
4/10
5/10
6/10
7/10
8/10
9/10
10/10

亲,该文档总共17页,到这已经超出免费预览范围,如果喜欢就直接下载吧~

如果您无法下载资料,请参考说明:

1、部分资料下载需要金币,请确保您的账户上有足够的金币

2、已购买过的文档,再次下载不重复扣费

3、资料包下载后请先用软件解压,在使用对应软件打开

实验报告 课程名称数值分析实验项目名称解线性方程组实验类型上机实验学时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.追赶法: 设系数矩阵为三对角矩阵 则方程组Ax=f称为三对角方程组。 设矩阵A非奇异,A有Crout分解A=LU,其中L为下三角矩阵,U为单位上三角矩阵,记 可先依次求出L,U中的元素后,令Ux=y,先求解下三角方程组Ly=f得出y,再求解上三角方程组Ux=y。 4.雅克比迭代法: 首先将方程组中的系数矩阵A分解成三部分,即:A=L+D+U,如图1所示,其中D为对角阵,L为下三角矩阵,U为上三角矩阵。 之后确定迭代格式,X=BX+f,如图2所示,其中B称为迭代矩阵,雅克比迭代法中一般记为J。(k=0,1,......)再选取初始迭代向量X,开始逐次迭代。 5.超松弛迭代法(SOR) 它是在GS法基础上为提高收敛速度,采用加权平均而得到的新算法。 选取分裂矩阵M为带参数的下三角矩阵 M=(D-), 其中>0为可选择的松弛因子,一般当1<<2时称为超松弛。 二.实验题目及实验目的 1.(第五章习题8)用直接三角分解(杜利特尔(Doolittle)分解)求线性方程组 ++=9, ++=8, ++=8 的解。 2.(第五章习题9)用追赶法解三对角方程组Ax=b,其中 A=,b=. 3.(第五章习题10)用改进的平方根法解线性方程组 = 4.(第六章习题7)用SOR方法解线性方程组(分别取松弛因子ω=1.03,ω=1,ω=1.1) 4- =1, -+4-=4, -+4=-3. 精确解x=(,1,-).要求当<5×10时迭代终止,并且对每一个ω值确定迭代次数. 5.(第六章习题8)用SOR方法解线性方程组(取ω=0.9) 5-2+=-12, -+4-2=20, 2-3+10=3. 要求当<10时迭代终止. 6.(第六章习题9)设有线性方程组Ax=b,其中A为对称正定阵,迭代公式 +ω(b-A),k=0,1,2…, 试证明当0<ω<时上述迭代法收敛(其中0<(A)). 7.(第六章计算实习题1)给出线性方程组Hx=b,其中系数矩阵H为希尔伯特矩阵: Hx=(h)R,h=,i,j=1,2,…,n. 假设x=(1,1,…,1)R,b=Hx.若取n=6,8,10,分别雅克比迭代法及SOR迭代(ω=1,1.25,1.5)求解.比较计算结果. 三.实验手段: 指操作环境和平台:win7系统下MATLABR2009a 程序语言:一种类似C语言的程序语言,但比C语言要宽松得多,非常方便。 四.程序 1. =1\*GB3①直接三角分解(文件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:1 h=s(i,m+1); forj=m:-1:1; ifj~=i h=h-x(j)*s(i,j); end end x(i)=h/s(i,i); end =2\*GB3②控制台输入代码: >>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. =1\*GB3①追赶法(文件ZG_SDJ.m) functionx=ZG_SDJ(a,b,c,f) %aÊǶԽÇÏßÔªËØ %bÊǶԽÇÏßÉÏ·½µÄÔªËØ£¬¸öÊý±ÈaÉÙÒ»¸ö %cÊǶԽÇÏßÏ·½µÄÔªËØ£¬¸öÊý±ÈaÉÙÒ»¸ö %fÊdz£ÊýÏîb N=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:N a1(j1)=a(j1)-c(j1)*b1(