数学实验 5：线性代数方程组的数值解法
线性代数方程组的数值解法讨论
线性代数方程组的数值解法讨论解线性方程组的方法，主要分为直接方法和迭代方法两种。
直接法是在没有舍入误差的假设下能在预定的运算次数内求得精确解。
而实际上，原始数据的误差和运算的舍入误差是不可以避免的，实际上获得的也是近似解。
迭代法是构造一定的递推格式，产生逼近精确解的序列。
对于高阶方程组，如一些偏微分方程数值求解中出现的方程组，采用直接法计算代价比较高，迭代法则简单又实用，因此比较受工程人员青睐。
小组成员本着工程应用，讨论将学习的理论知识转变为matlab 代码。
讨论的成果也以各种代码的形式在下面展现。
1 Jacobi 迭代法使用Jacobi 迭代法，首先必须给定初始值，其计算过程可以用以下步骤描述： 步骤1 输入系数矩阵A ，常熟向量b ，初值(0)x ，误差限ε，正整数N ，令1k =.步骤2 (0)11ni i ij jj ii j i x b a x a =≠⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦∑，(0)j x 代表(0)x 的第j 个分量。
步骤3 计算11ni i ij j j ii j i y b a x a =≠⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦∑，判断1max i i i n x y ε≤≤-<，如果是，则结束迭代，转入步骤5；否则，转入步骤4。
步骤4 判断k N =？如果是，则输出失败标志；否则，置1k k =+，i i x y ⇐，1,2,,i n =，转入步骤2。
步骤5 输出12,,n y y y 。
雅可比迭代代码function [x,k]=Fjacobi(A,b,x0,tol)% jacobi 迭代法 计算线性方程组% tol 为输入误差容限，x0为迭代初值max1= 300; %默认最多迭代300，超过要300次给出警告 D=diag(diag(A)); L=-tril(A,-1);U=-triu(A,1); B=D\(L+U); f=D\b; x=B*x0+f;k=1; %迭代次数while norm(x-x0)>=tol x0=x;x=B*x0+f; k=k+1;if(k>=max1)disp('迭代超过300次，方程组可能不收敛'); return; end%[k x'] %显示每一步迭代的结果 End2 高斯赛德尔迭代由Jacobi 迭代法中，每一次的迭代只用到前一次的迭代值，若每一次迭代充分利用当前最新的迭代值，即在计算第i 个分量(1)k i x +时，用最新分量11()k x +，12()k x +…(1)1k i x +-代替旧分量)1(k x ', )2(k x …)3(k x 就得到高斯赛德尔迭代格式，其数学表达式为：1(1)(1)()111(1,2,,)i n k k k ii ij j ij j j j i ii xb a x a x i n a -++==+⎛⎫=--= ⎪⎝⎭∑∑具体形式如下：()()()(1)()()()11221331111(1)(1)()()22112332222(1)(1)(1)(1)(1)112233,11111k k k k n n k k k k n n k k k k k n n n n n n n n nnx a x a x a x b a x a x a x a x b a x a x a x a x a x b a ++++++++--=----+=----+⋯⋯⋯⋯⋯⋯=-----+矩阵形式表示为：()(1)1(1)()(0,1,2,,),k k k k n +-+=++=x D Lx Ux b将(1)(1)()(0,1,2,,)k k k k n ++=++=Dx Lx Ux b 移项整理得： (1)1()1()()(0,1,2,,))k k x D L Ux D L b k n +--=-+-=记11(),()--=-=-M D L U g D L b ，则(1)()k k x x +=+M g高斯塞德尔迭代function [x,k]=Fgseid(A,b,x0,tol)%高斯-塞德尔迭代法 计算线性方程组 % tol 为误差容限max1= 300; %默认最高迭代300次D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); G=(D-L)\U; f=(D-L)\b; x=G*x0+f;k=1; while norm(x-x0)>=tol x0=x;x=G*x0+f; k=k+1;if(k>=max1)disp('迭代次数太多，可能不收敛'); return; end% [k,x'] %显示每一步迭代结果 End3 超松弛迭代法在工程中最常遇到的问题便是线性代数方程组的求解,而线性代数方程组的求解一般可以分为两类,一类是直接法(精确法),包括克莱姆法则方法、LD 分解法等，另一类是迭代法（近似法），包括雅克比迭代法、高斯迭代法、超松弛迭代法等。
实验5_线性代数方程组的数值解法
实验5 线性代数方程组的数值解法化工系 毕啸天 2010011811【实验目的】1. 学会用MATLAB 软件数值求解线性代数方程组，对迭代法的收敛性和解的稳定性作初步分析；2. 通过实例学习用线性代数方程组解决简化的实际问题。
【实验内容】题目3已知方程组Ax=b ，其中A ，定义为⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--------------=32/14/12/132/14/14/12/132/14/14/12/132/14/12/13A试通过迭代法求解此方程组，认识迭代法收敛的含义以及迭代初值和方程组系数矩阵性质对收敛速度的影响。
实验要求：（1）选取不同的初始向量x (0)和不同的方程组右端项向量b ，给定迭代误差要求，用雅可比迭代法和高斯-赛德尔迭代法计算，观测得到的迭代向量序列是否均收敛？若收敛，记录迭代次数，分析计算结果并得出你的结论；（2）取定右端向量b 和初始向量x (0)，将A 的主对角线元素成倍增长若干次，非主对角线元素不变，每次用雅可比迭代法计算，要求迭代误差满足 ，比较收敛速度，分析现象并得出你的结论。
3.1 模型分析选取初始向量x(0) =(1,1,…,1)T ，b=(1,1,…,1)T ，迭代要求为误差满足 ，编写雅各比、高斯-赛德尔迭代法的函数，迭代求解。
3.2 程序代码function x = Jacobi( x0,A,b,m ) D=diag(diag(A)); U=-triu(A,1); L=-tril(A,-1); B1=D\(L+U); f1=D\b; x(:,1)=x0;x(:,2)=B1*x(:,1)+f1; k=1;while norm((x(:,k+1)-x(:,k)),inf)>m x(:,k+2)=B1*x(:,k+1)+f1; k=k+1; endendfunction x = Gauss( x0,A,b,m )D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B2=(D-L)\U;f2=(D-L)\b;x(:,1)=x0;x(:,2)=B2*x(:,1)+f2;k=1;while norm((x(:,k+1)-x(:,k)),inf)>mx(:,k+2)=B2*x(:,k+1)+f2;k=k+1;endendA1=3.*eye(20,20);A2=sparse(1:19,2:20,-1/2,20,20);A3=sparse(1:18,3:20,-1/4,20,20);AA=A1+A2+A3+A2'+A3';A=full(AA);b=ones(20,1); %输入自选右端项向量b x0=ones(20,1); %输入自选初始向量x0 m=1e-5;x1=Jacob(x0,A,b,m);x2=Gauss(x0,A,b,m);结果输出数据：3.3.1将x0各分量初值置为0【分析】从数据中可以看出，当迭代的初值变化了，达到相同精度所需要的迭代次数也变化了。
线性方程组的数值解法实验报告
实验报告——线性方程组的数值解法姓名：班级：学号：日期：一 实践目的1. 熟悉求解线性方程组的有关理论和方法。
2. 会编列主元消去法，全主元消去法，雅克比迭代法及高斯-赛德尔迭代法的程序。
3.通过实际计算，进一步了解各种方法的优缺点，选择合适的数值方法。
4. 进一步应用数学知识，扩展数学思维，提高编程能力。
二 问题定义及题目分析1.求解线性方程是实际中常遇到的问题。
而一般求解方法有两种，一个是直接求解法，一个是迭代法。
2.直接法是就是通过有限步四则运算求的方程准确解的方法。
但实际计算中必然存在舍入误差，因此这种方法只能得到近似解。
3.迭代法是先给一个解的初始近似值，然后按一定的法则求出更准确的解，即是用某种极限过程逐步逼近准确解的方法。
4.这次实践，将采用直接法：列主元高斯消去法，全主元高斯消去法；迭代法：雅克比迭代法和高斯-赛德尔迭代法。
三 详细设计 设有线性方程组11112211n n a x a x a x b +++=21122222n n a x a x a x b +++= 1122n n nn n n a x a x a x b +++=令A= 111212122212n n n n nn a a a aa a a a a ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦ x=12n x x x ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦ b= 12n b b b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦一. 上三角方程组的求解很简单。
所以若能将方程组化为上三角方程组，那就很容易得到方程的解，高斯消去法则是一种简洁高效的方法，可以将方程组化为上三角方程组，但是这种方法必须满足每一步时0kk a =才能求解，还有当kk a 绝对值很小时，作分母会引起较大的舍入误差。
因此在消元过程中应尽量选择绝对值较大的系数作为主元素。
1.列主元消去法很好的解决这一问题。
将1x 的系数1(1)k a k n ≤≤中绝对值最大者作为主元素，交换第一行和此元素所在的行，然后主元消去非主元项1x 的系数，。
实验五(线性方程组的数值解法和非线性方程求解)
1大学数学实验 实验报告 | 2014/4/5一、 实验目的1、学习用Matlab 软件数值求解线性代数方程组，对迭代法的收敛性和解的稳定性作初步分析；2、通过实例学习用线性代数方程组解决简化问题。
二、 实验内容项目一：种群的繁殖与稳定收获：种群的数量因繁殖而增加，因自然死亡而减少，对于人工饲养的种群（比如家畜）而言，为了保证稳定的收获，各个年龄的种群数量应维持不变。
种群因雌性个体的繁殖而改变，为方便起见以下种群数量均指其中的雌性。
种群年龄记作k=1,2,…,n ,当年年龄k 的种群数量记作x k ，繁殖率记作b k （每个雌性个体1年的繁殖的数量），自然存活率记作s k （s k =1−d k ，d k 为1年的死亡率），收获量记作ℎk ，则来年年龄k 的种群数量x ̌k 应该为x ̌k =∑b k n k=1x k , x ̌k+1=s k x k −ℎk , (k=1,2,…,n -1)。
要求各个年龄的种群数量每年维持不变就是要求使得x ̌k =x k , (k=1,2,…,n -1).(1) 如果b k , s k 已知，给定收获量ℎk ，建立求各个年龄的稳定种群数量x k 的模型（用矩阵、向量表示）.(2) 设n =5,b 1=b 2=b 5=0,b 3=5,b 4=3,s 1=s 4=0.4,s 2=s 3=0.6,如要求ℎ1~ℎ5为500,400,200,100,100，求x 1~x 5.(3) 要使ℎ1~ℎ5均为500，如何达到？问题分析：该问题属于简单的种群数量增长模型，在一定的条件（存活率，繁殖率等）下为使各年龄阶段的种群数量保持不变，各个年龄段的种群数量将会满足一定的要求，只要找到种群数量与各个参量之间的关系，建立起种群数量恒定的方程就可以求解出各年龄阶段的种群数量。
模型建立：根据题目中的信息，令x ̌k =x k ，得到方程组如下：{x ̌1=∑b k nk=1x k =x 1x ̌k+1=s k x k −ℎk =x k+1整理得到：{−x 1∑b k nk=1x k =0−x k+1+s k x k =ℎk2 大学数学实验 实验报告 | 2014/4/52写成系数矩阵的形式如下：A =[b 1−1b 2b 3s 1−100s 2−1…b n−1b n0000⋮⋱⋮000000000⋯00−10s n−1−1]令h =[0, ℎ1,ℎ2,ℎ3,…,ℎn−2,ℎn−1]Tx =[x n , x n−1,…,x 1]T则方程组化为矩阵形式：Ax =h ，即为所求模型。
线性方程组的数值解法
对每行计算乘数
mi1aa1i1111, i2,3,,n
用 mi1 乘以第1个方程，加到第 i个方程，消去
第 2个方程到第 n个方程的未知数x1 ，得 A2xb2
即:
a111
a112 a222
aa1212nnxx12
bb1212
an22 an2nxn bn2
其中： a bii2 2 j a bii1 1j m m ii1 1b a1 1 1 1j i,j2,3, ,n
an1 ann
x1
x
x n
b1
b
b n
若矩阵A非奇异，即A的行列式 deAt0，根据
克莱姆（Gramer）法则，方程组有唯一 解：
xi
Di D
i1,2, ,n
其中D表示 detA，D i 表示 D 中第 i列换成 b后
所得的行列式。
当阶数较高时用这种方法求解是不现实的。n阶行
综上所述，高斯消去法的框图如图3-1所示。从 中可看出高斯消去法的计算机运算和存储方式的特点：
1〉按消元规则进行运算后，对角线以下元素为0。 故对于对角线以下元素不用作计算，减小了计算量。
2〉对角线以下元素对回代求解无影响，故可将乘 数放在该处，即
a akikkaik,ik1,k2, ,n
以节省存储单元。
列式有 n项!，每项又是 个n数的乘积。对较大的 ，
其计n算量之大，是一般计算机难以完成的。而且， 这时的舍入误差对计算结果的影响也较大。
例如，求解一个20阶线性方程组，用加减消元法需 3000次乘法运算，而用克莱姆法则要进行 9.71020次 运算，如用每秒1亿次乘法运算的计算机要30万年。
线性代数方程组的计算机解法常用方法：
线性代数方程组的数值解法_百度文库
线性代数方程组的数值解法【实验目的】1. 学会用MATLAB 软件数值求解线性代数方程组，对迭代法的收敛性和解的稳定性作初步分析；2. 通过实例学习用线性代数方程组解决简化的实际问题。
【实验内容】【题目1】通过求解线性方程组A1x=b1和A2x=b2，理解条件数的意义和方程组的性态对解的影响。
其中A1是n阶范德蒙矩阵，即⎡1x0⎢1x1⎢A1=⎢⎢⎢⎣1xn-12x0x12 2xn-1n-1⎤ x0⎥ x1n-1⎥1,...,n-1 ，xk=1+0.1k，k=0,⎥ n-1⎥ xn-1⎥⎦A2是n阶希尔伯特矩阵，b1，b2分别是A1，A2的行和。
（1）编程构造A1（A2可直接用命令产生）和b1，b2；你能预先知道方程组A1x=和A2x=。
b2的解吗？令n=5，用左除命令求解（用预先知道的解可检验程序）b1（2）令n=5,7,9,…，计算A1，A2的条件数。
为观察它们是否病态，做以下试验：b1，b2不变，A1和A2的元素A1(n,n)，A2(n,n)分别加扰动ε后求解；A1和A2不变，b1，b2的分量b1(n)，分析A和b的微小扰动对解的影响。
b2(n)分别加扰动ε求解。
ε取10-1010，-8，10-6。
（3）经扰动得到的解记做x~，计算误差-x~x，与用条件数估计的误差相比较。
1.1构造A1，A2和b1，b2首先令n=5，构造出A1，A2和b1，b2。
首先运行以下程序，输出A1。
运行以下程序对A1,A2求行和：由于b1，b2分别是A1，A2的行和，所以可以预知x1=运行下列程序，用左除命令对b1，b2进行求解：得到以下结果： T。
x2=(1,1, ,1)1.2 计算条件数并观察是否为病态1.不加扰动，计算条件数。
运行以下程序：由此可知，A1，A2的条件数分别是3.574∗10, 4,766∗10。
2.b1，b2不变，A1(n,n)，A2(n,n)分别加扰动（1）n=5时设x11,x12,x13分别为A1添加扰动10−10，10−8,10−6后的解。
数值分析实验报告-清华大学--线性代数方程组的数值解法
数值分析实验报告-清华大学--线性代数方程组的数值解法(总15页)--本页仅作为文档封面，使用时请直接删除即可----内页可以根据需求调整合适字体及大小--线性代数方程组的数值解法实验1. 主元的选取与算法的稳定性问题提出：Gauss 消去法是我们在线性代数中已经熟悉的。
但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss 消去法作为数值算法的稳定性呢？Gauss 消去法从理论算法到数值算法，其关键是主元的选择。
主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。
实验内容：考虑线性方程组 n n n R b R A b Ax ∈∈=⨯,,编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss 消去过程。
实验要求：（1）取矩阵⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=1415157,6816816816 b A ，则方程有解T x )1,,1,1(* =。
取n=10计算矩阵的条件数。
让程序自动选取主元，结果如何？（2）现选择程序中手动选取主元的功能。
每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。
若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。
（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。
（4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。
重复上述实验，观察记录并分析实验结果。
程序清单n=input('矩阵A 的阶数:n=');A=6*diag(ones(1,n))+diag(ones(1,n-1),1)+8*diag(ones(1,n-1),-1); b=A*ones(n,1);p=input('计算条件数使用p-范数,p='); cond_A=cond(A,p) [m,n]=size(A);Ab=[A b];r=input('选主元方式(0:自动;1:手动)，r=');Abfor i=1:n-1switch rcase(0)[aii,ip]=max(abs(Ab(i:n,i)));ip=ip+i-1;case (1)ip=input(['第',num2str(i),'步消元，请输入第',num2str(i),'列所选元素所处的行数：']);end;Ab([i ip],:)=Ab([ip i],:);aii=Ab(i,i);for k=i+1:nAb(k,i:n+1)=Ab(k,i:n+1)-(Ab(k,i)/aii)*Ab(i,i:n+1);end;if r==1Abendend;x=zeros(n,1);x(n)=Ab(n,n+1)/Ab(n,n);for i=n-1:-1:1x(i)=(Ab(i,n+1)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);endx运行结果(1)n=10，矩阵的条件数及自动选主元Cond(A,1) =×103Cond(A,2) = ×103Cond(A,inf) =×103程序自动选择主元（列主元）a.输入数据矩阵A的阶数:n=10计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动)，r=0b.计算结果x=[1，1，1，1，1，1，1，1，1，1]T(2)n=10，手动选主元a. 每步消去过程总选取按模最小或按模尽可能小的元素作为主元矩阵A 的阶数:n=10计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动)，r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元，请输入第1列所选元素所处的行数：1(2)(2) 6.0000 1.00007.00004.6667 1.0000 5.66678.0000 6.000015.0000[]8.00001.000015.00006.0000 1.00008.0000 6.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元，请输入第2列所选元素所处的行数：2…（实际选择时，第k 步选择主元处于第k 行） 最终计算得x=[, , , , , , , , , ]Tb. 每步消去过程总选取按模最大的元素作为主元 矩阵A 的阶数:n=10计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动)，r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元，请输入第1列所选元素所处的行数：2(2)(2)8.0000 6.0000 1.000015.0000-3.50000.7500-4.250008.0000 6.0000 1.000015.0000[]8.0000 6.000015.00008.0000 1.00006.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元，请输入第2列所选元素所处的行数：3…（实际选择时，第k 步选择主元处于第k+1行） 最终计算得x=[1，1，1，1，1，1，1，1，1，1]T(3)n=20，手动选主元a. 每步消去过程总选取按模最小或按模尽可能小的元素作为主元 矩阵A 的阶数:n=20计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动)，r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元，请输入第1列所选元素所处的行数：1(2)(2) 6.0000 1.00007.00004.6667 1.0000 5.66678.0000 6.000015.0000[]8.00001.000015.00006.0000 1.00008.0000 6.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元，请输入第2列所选元素所处的行数：2…（实际选择时，第k 步选择主元处于第k 行） 最终计算得x=[，，，，，，，，，，，，，，，，，，，]T b. 每步消去过程总选取按模最大的元素作为主元 矩阵A 的阶数:n=20计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动)，r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元，请输入第1列所选元素所处的行数：2(2)(2)8.0000 6.0000 1.000015.0000-3.50000.7500-4.250008.0000 6.0000 1.000015.0000[]8.0000 6.000015.00008.0000 1.00006.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元，请输入第2列所选元素所处的行数：3…（实际选择时，第k步选择主元处于第k+1行）最终计算得x=[1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1]T(4)A分别为幻方矩阵，Hilbert矩阵，pascal矩阵和随机矩阵简要分析计算(1)表明：对于同一矩阵，不同范数定义的条件数是不同的；Gauss消去法在消去过程中选择模最大的主元能够得到比较精确的解。
计算方法实验报告-线性方程组的数值解法
重庆大学学生实验报告实验课程名称计算方法开课实验室DS1421学院年级专业学生姓名学号开课时间至学年第学期1．实验目的（1）高斯列主元消去法求解线性方程组的过程（2）熟悉用迭代法求解线性方程组的过程（3）设计出相应的算法，编制相应的函数子程序2．实验内容分别用高斯列主元消去法 ,Jacobi 迭代法,Gauss--Saidel 迭代法,超松弛迭代法求解线性方程组⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡------725101391444321131243301024321x x x x 3.实验过程解：（1）高斯列主元消去法编制高斯列主元消去法的M 文件程序如下：%高斯列主元消元法求解线性方程组Ax=b%A 为输入矩阵系数,b 为方程组右端系数%方程组的解保存在x 变量中format long;%设置为长格式显示,显示15位小数A=[2,10,0,-3;-3,-4,-12,13;1,2,3,-4;4,14,9,-13]b=[10,5,-2,7]'[m,n]=size(A);%先检查系数正确性if m~=nerror('矩阵A 的行数和列数必须相同');return;endif m~=size(b)error('b 的大小必须和A 的行数或A 的列数相同');return;end%再检查方程是否存在唯一解if rank(A)~=rank([A,b])error('A 矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解');return;endc=n+1;A(:,c)=b; %(增广)for k=1:n-1[r,m]=max(abs(A(k:n,k))); %选主元m=m+k-1; %修正操作行的值if(A(m,k)~=0)if(m~=k)A([k m],:)=A([m k],:); %换行endA(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k))*A(k, k:c); %消去 endendx=zeros(length(b),1); %回代求解x(n)=A(n,c)/A(n,n);for k=n-1:-1:1x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n))/A(k,k);enddisp('X=');disp(x);format short;%设置为默认格式显示,显示5位运行，结果如下所示：（2）Jacobi迭代法编制迭代计算的M文件程序如下:%Jacobi迭代法求解% A为方程组的增广矩阵clc;A=[2,10,0,-3,10;-3,-4,-12,13,5;1,2,3,-4,-2;4,14,9,-13,7] MAXTIME=50;%最多进行50次迭代eps=1e-5;%迭代误差[n,m]=size(A);x=zeros(n,1);%迭代初值y=zeros(n,1);k=0;%进入迭代计算disp('迭代过程X的值情况如下:')disp('X=');while 1disp(x');for i=1:1:ns=0.0;for j=1:1:nif j~=is=s+A(i,j)*x(j);endy(i)=(A(i,n+1)-s)/A(i,i);endendfor i=1:1:nmaxeps=max(0,abs(x(i)-y(i))); %检查是否满足迭代精度要求 endif maxeps<=eps%小于迭代精度退出迭代for i=1:1:nx(i)=y(i);%将结果赋给xendreturn;endfor i=1:1:n%若不满足迭代精度要求继续进行迭代x(i)=y(i);y(i)=0.0;endk=k+1;if k>MAXTIME%超过最大迭代次数退出error('超过最大迭代次数,退出');return;endend运行该程序结果如下:（3）Gauss--Saidel迭代法编制求解程序Gauss_Seidel.m如下:%Gauss_Seidel.m% A为方程组的增广矩阵clc;format long;A=[2,10,0,-3,10;-3,-4,-12,13,5;1,2,3,-4,-2;4,14,9,-13,7][n,m]=size(A);%最多进行50次迭代Maxtime=50;%控制误差Eps=10E-5;%初始迭代值x=zeros(1,n);disp('x=');%迭代次数小于最大迭代次数,进入迭代for k=1:Maxtimedisp(x);for i=1:ns=0.0;for j=1:nif i~=js=s+A(i,j)*x(j);%计算和endendx(i)=(A(i,n+1)-s)/A(i,i);%求出此时迭代的值end%因为方程的精确解为整数,所以这里将迭代结果向整数靠近的误差作为判断迭代是否停止的条件if sum((x-floor(x)).^2)<Epsbreak;end;end;X=x;disp('迭代结果:');Xformat short;运行结果如下所示：（4）超松弛迭代法编写函数M文件如下：%SOR法求解%w=1.45%方程组系数矩阵clc;A=[2,10,0,-3;-3,-4,-12,13;1,2,3,-4;4,14,9,-13]%方程组右端系数b=[10,5,-2,7]'w=1.45;%最大迭代次数Maxtime=100;%精度要求Eps=1E-5;%以15位小数显示format long;n=length(A);k=0;%初始迭代值x=ones(n,1);y=x;disp('迭代过程:');disp('x=');while 1y=x;disp(x');%计算过程for i=1:ns=b(i);for j=1:nif j~=is=s-A(i,j)*x(j);endendif abs(A(i,i))<1E-10 | k>=Maxtimeerror('已达最大迭代次数或矩阵系数近似为0,无法进行迭代'); return;ends=s/A(i,i);x(i)=(1-w)*x(i)+w*s;endif norm(y-x,inf)<Eps%达到精度要求退出计算break;endk=k+1;enddisp('最后迭代结果:');%最后的结果X=x'%设为默认显示格式format short; 结果如下：4．实验环境及实验文件存档名实验环境：Matlab7.0文件存档名：Gauss.m，Jacobi.m，Gauss_Seidel.m，SOR.m5.实验结果及分析=1.0000=2.0000=3.0000=4.0000经过验证，高斯列主元消结果正确。
线性方程组的数值解法-安振华-2012011837
实验5：线性方程组的数值解法化学工程系分2 安振华2012011837【实验目的】1、掌握线性方程组的常用数值解法，包括高斯消去法、LU分解法以及校正法。
2、体验数值计算的时间复杂度和计算规模的关系。
3、加深对数值计算误差的理解。
4、学习使用迭代法等算法，求解非线性方程。
5、学习如何使用MATLAB解非线性方程组和方程组。
【实验容】【实验五：习题9】种群的繁殖与稳定收获：种群的数量因繁殖而增加，因自然死亡而减少，对于人工饲养的种群（比如家畜）而言，为了保证稳定的收获，各个年龄的种群数量应保持不变，种群因雌性个体的繁殖而改变，为方便起见以下种群数量均指其中的雌性。
种群年龄记作k=1,2,…,n，当年年龄k的种群数量记作x k，繁殖率记作b k（每个雌性个体在1年繁殖的数量），自然存活率记作s k（s k=1-d k，d k为1年的死亡率），收获量记作h k，则来年年龄k的种群数量k x应为：111,(1,2,,1)n k k k k k k k x b x x s x h k n +===-=⋅⋅⋅-∑要求各个年龄的种群数量每年维持不变就是要使(1,2,,)k k x x k n ==⋅⋅⋅（1） 若b k ，s k 已知，给定收获量h k ，建立求各年龄的稳定种群数量x k 的模型（用矩阵向量表示）（2） 设n=5，b 1=b 2=b 5=0，b 3=5，b 4=3，s 1=s 4=0.4，s 2=s 3=0.6，如果要求h 1~h 5为500,400,200,100,100,求x 1~x 5 （3） 要使h 1~h 5均为500，如何达到？ 【分析】为方便起见以下种群数量均指其中的雌性。
我们并且有以下的假设：(1)雌性个体的繁殖率和存活率在特定的时间是不变的。
(2)人工饲养的种群在质量和数量上是不受外界环境和资源的限制的。
(3)模型中不考虑人为的或是自然的灾害所造成的种群数量、繁殖率和存活率的变动。
数学建模线性方程组的数值解法
直接法 经过有限次算术运算求出精确解（实际上 由 于 有 舍 入 误 差 只 能 得 到 近 似 解 ） ---- 高 斯 （Gauss）消元法及与它密切相关的矩阵LU分解 迭代法 从初始解出发，根据设计好的步骤用逐次 求出的近似解逼近精确解 ---- 雅可比（Jacobi） 迭代法和高斯—塞德尔（Gauss—Seidel）迭代法
(k )
0.1x1
( k 1)
0.3x2
( k 1)
1.4
Gauss-Seideil迭代公式 Dx ( k 1) Lx ( k 1) Ux ( k ) b
用它作除数会导致舍入误 差的很大增加 解决 办法 选
(k ) aik
(i k , n) 最大的一个(列主元)
将列主元所在行与第k行交换后, 再按上面的 高斯消元法进行下去，称为列主元消元法。
直接法 - 高斯消元法的矩阵表示
高斯消元法的第一次消元
a11 x1 a12 x2 a1n xn b1 a21 x1 a22 x2 a2 n xn b2 an1 x1 an 2 x2 ann xn bn
数值解法(迭代解法)的收敛性
实验5的主要内容
1. 两类数值解法: 直接方法；迭代方法
2. 超定线性方程组的最小二乘解 3. 线性方程组数值解法的MATLAB实现 4. 实际问题中方程组的数值解
线性方程组的一般形式、两类解法
a11 x1 a12 x 2 a1n x n b1 a 21 x1 a 22 x 2 a 2 n x n b2 a n1 x1 a n 2 x 2 a nn x n bn
大学数学实验
Mathematical Experiments 实验5 线性代数方程组的数值解法
数学实验 5：线性代数方程组的数值解法
大倍数 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 9.0000 10.0000 11.0000 12.0000 13.0000 14.0000 15.0000 16.0000 17.0000 18.0000 19.0000 20.0000
小值q 0.4893 0.2447 0.1631 0.1223 0.0979 0.0816 0.0699 0.0612 0.0544 0.0489 0.0445 0.0408 0.0376 0.0350 0.0326 0.0306 0.0288 0.0272 0.0258 0.0245 21.0000 12.0000 9.0000 8.0000 8.0000 7.0000 7.0000 7.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 5.0000 5.0000
实验 5：线性代数方程组的数值解法
习题3：
已知方程组,其中，定义为： 试通过迭代法求解此方程组，认识迭代法收敛的含义以及迭代初值和方 程组系数矩阵性质对收敛速度的影响。实验要求： （1） 选取不同的初始向量x0和不同的方程组右端向量b，给定迭 代误差要求，用雅可比迭代法和高斯-赛德尔迭代法计算， 观测得到的迭代向量序列是否均收敛？若收敛，记录迭代 次数，分析计算结果并得出结论； （2） 取定右端向量b和初始向量x0，将A的主对角线元素成倍的 增长若干次，非主对角元素不变，每次用雅可比迭代法计 算，要求迭代误差满足,比较收敛速度，分析现象并得出结 论。 1、 程序设计（可直接粘贴运行） 1） Jacobi迭代法 function y=jacobi(a,b,x0,e,m) %定义jacobi函数，其中：a,b为线性方程组中的矩阵和右端向量；x0 为初始值； %e和m分别为人为设定的精度和预计迭代次数；运行结果y为迭代的结 果和所有中间值组成的 %矩阵 y=0; %对y初始化 d=diag(diag(a)); %按雅可比迭代标准形 形式取主对角元素作为矩阵D u=-triu(a,1); %取上三角矩阵u l=-tril(a,-1); %取下三角矩阵l bj=d^-1*(l+u); fj=d^-1*b; x=[x0,zeros(20,m-1)]; %初始化x,其中x1=x0,即 初始值 for k=1:m %人为规定迭代次 数，防止不收敛迭代导致死循环 x(:,k+1)=bj*x(:,k)+fj; %jacobi迭代 if norm(x(:,k+1)-x(:,k),inf)<e
应用数值分析线性代数方程组数值解法
线性代数方程组数值解法——直接法一、n 阶线性代数方程组b x A =⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡n n nn n n n n b b b x x x a a a a a a a a a 2121212222111211二、上（下）方程组与回代（前推）过程 1、上三角方程组⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡n n nn n n b b b x x x a a a a a a 2121222112112、回代过程⎪⎩⎪⎨⎧-=⎪⎭⎫ ⎝⎛-==∑+=)1,2,1(,1 n i a x a b x a b x ii ni j j ij ii nn n n 3、下三角方程组⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡n n nn n n b b b x x x a a a a a a 2121212221114、前推过程⎪⎩⎪⎨⎧=⎪⎭⎫ ⎝⎛-==∑-=),2,1(,111111n i a x a b x a b x ii i j j ij ii 三、顺序Gauss 消去法 1、消元过程（1）n 阶方程组经1-k 次消元后，得到的新的系数矩阵)(k A 具有如下形式：⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡=)()()()()2(2)2(22)1(1)1(12)1(11)(k nn k nk k knk kk n n k a a a a a a a a a A（2）在第k 次，计算乘数),1(,)()(n k i a a m k kk k ik ik +==，则)1(+k A元素的计算公式为：⎩⎨⎧+=-=+=-=++),1(,),1,(,)()()1()()()1(n k i b m b b n k j i a m a a k k ik k i k ik kj ik k ij k ij （3）完成1-n 次消元后：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡)()2(2)1(121)()2(2)2(22)1(1112111n n n n nn n nb b b x x x a a a a a a ）（）（ 2、回代过程⎪⎩⎪⎨⎧-=⎪⎭⎫ ⎝⎛-==∑+=)1,2,1(,)(1)()()()( n k a x a b x a b x k kk n k j j k kj k k k n nn n n n3、计算量估计（1）第k 次消元需要进行)1)((+--k n k n 次乘法和)(k n -次除法：6523)1()()(231111nn n k n k n k n n k n k -+=+--+-∑∑-=-= （2）回代过程乘除法次数：22)(211nn n k n n k +=+-∑-= （3）乘数法总次数：3323n n n MD -+=（4）加减法总次数：652323n n n AS -+=4、要求：系数矩阵A 的所有顺序主子式),1,2,1(0n n k k -=≠∆ 四、列主元Gauss 消去法1、对于非奇异矩阵A ，在第一步消元时，即使A 的第一个元素为零，但是A 的第一列元素中至少有一个不为零，把这个不为零的元素所在的行与第一行交换位置，即可进行消元。
数学实验第3次作业_线性方程组的数值解法
线性方程组的数值解法一 实验目的1 学会用MATLAB 软件求解线性代数方程组，对迭代法的收敛性和解的稳定性做初步分析；2 通过实例学习线性代数方程组解决简化的实际问题。
二 实验内容1 已知方程组Ax =b ，其中A ∈ℝ20×20，定义为：31/21/41/231/21/41/41/231/21/41/41/231/21/41/23--⎡⎤⎢⎥---⎢⎥⎢⎥---⎢⎥-⎢⎥⎢⎥---⎢⎥--⎣⎦试通过迭代法求解此方程组，认识迭代法收敛的含义以及迭代初值和方程组系数矩阵性质对收敛速度的影响。
实验要求：(1) 选取不同的初始向量x (0)和不同的方程组右端项向量b ，给定迭代误差要求，用雅可比迭代法和高斯-赛德尔迭代法计算，观测得到的迭代向量序列是否均收敛？若收敛，记录迭代次数，分析计算结果并得出你的结论。
(2) 取定右端项向量b 和初始向量x (0)，将A 的主对角线元素成倍增长若干次，非对角线元素不变，每次用雅可比迭代法计算，要求迭代误差满足‖x (k+1)−x (k)‖∞<10−5，比较收敛速度，分析现象并得出你的结论。
初步解决：首先建立利用雅克比迭代法和高斯-赛德尔迭代法计算的函数M 文件。
雅克比迭代法：高斯-赛德尔迭代法：关于迭代法的收敛性：由原式Ax=b化简可得x=Bx+f，也就是说，研究此处迭代法是否收敛实际上就是研究有原矩阵A经变换之后所得矩阵B的性质。
而由于该题第一问中，A保持不变，所以对于收敛性只需研究一次。
在命令栏中输入以下命令：(1)构造题目中的矩阵A：(2)提取对角矩阵、上三角矩阵、下三角矩阵：(3)对于两种迭代方法分别求B：(4)分别计算两个矩阵的三种范数：n1和m1、n2和m2、n3和m3分别是矩阵B1、B2的1-范数、2-范数和∞范数。
由书上定理，矩阵的谱半径不超过任何一个范数，即ρ(B)≤‖B‖，而由图中可以看出，两个矩阵的六个范数没有一个大于1，所以两个矩阵的谱半径一定都小于1，所以此时两种迭代方法均收敛。
【清华】7.0_实验5-线性代数方程组的数值解法
0
0.5000
0.7500
0.9926
0.99979 0.99999 0.99999
0
0.5000
0.7500
0.9930
0.99981 0.99999 0.99999
0
0.5000
0.7500
0.9935
0.99983 0.99999
1
0
0.5000
0.7569
0.9943
0.99985 0.99999
B1=D\(L+U);
f1=D\b; x(:,1)=x0; x(:,2)=B1*x(:,1)+f1; k=1; while norm((x(:,k+1)-x(:,k)),inf)>m
x(:,k+2)=B1*x(:,k+1)+f1; k=k+1; end
%赋初值 %以差值向量的行范数作为误差判断标准
0.22222
0 0.083333 0.095238
0.21996
0 0.083333 0.095016
0 0.33333
0.6521
0 0.16667
0.22152
0 0.083333 0.095181
0 0.33333
0.66094
0 0.16667
0.22207
0 0.083333 0.095231
0 0.33333
0.6643
0 0.16667
0
0.7083
0.9069
0.9967
0.99999
1
1
0
0.6806
0.8960
0.9962
0.99998 0.99999
线性代数方程组的数值解法
) a
a / a (k1) kj
( k 1) kk
a a (k1)
ij
( k 1) ik
( j) kj
( j k 1,, n 1) (i k 1,, n; j k 1,n 1)
第n步：得到：
1
a (1) 12
a (1) 13
a (1) 1n
1
a (2) 23
a (2) 2n
a (1) 1n 1
A b 行初等变换 I x
相应地，计算公式可表述为：
对 k 1,2,, n, 依次计算
aaik((jjkk
) )
a(k 1) kj
a(k 1) ij
/
a(k kk
1)
a a (k 1) (k )
ik
kj
(
j
k
1, k
2,, n 1)
(i 1,2,, k 1, k 1, k 2,, n; j k 1,, n 1)
后用第i行元数（i
2,
,
n)减去第一行对应元素的a
(0) i1
倍，(i 2,, n),这样，a1(01)位置变为1，其余各行的第一
个元素变为0。
1
(A(0)
,
b(0)
)
a
(0) 21
a
(0) n1
a(0) 12
/
a(0) 11
a(0) 22
a(0) 13
/
a(0) 11
a(0) 23
a(0) 1n
0
a (1) 13
a (2) 23
a (1) 1n 1
a (2) 2n 1
记为
a (2) 33
a
(2) 3n 1
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
11.0000 5.8056 3.7199 2.7801 2.3626 2.1717 … 2.0072 2.0072
Gauss K=1 X1 X2 X3 5.3333 6.5556 7.5370
2.0540 1.1354 0.8539 0.7657 0.7377 … 0.7247 0.7247 2.9432 1.8459 1.5033 1.3949 1.3605 … 1.3444 1.3444 3.7386 2.5553 2.1815 2.0627 2.0249 … 2.0072 2.0072
结果分析： 根据判断迭代是否收敛的原理，若A是严格对角占优的，即，则 Jacobi迭代和Gauss迭代均收敛。由此推测，A的对角占优程度可能是影 响收敛速度的关键因素。由上述试验表面，对A的对角元素数值的扩大 的确可以明显的改善Jacobi迭代的收敛速度，与预测的情况相符。 主对角占优程度可以影响收敛速度，但是这个结论因为数学水平有
2、 运行结果分析 1） 不同初值（x0）和右端项（b）条件下的解的情况
表一：收敛性判断
1-范数 Jacobi 0.0163 0.0167 2-范数 0.0163 -范数 0.0167
Gauss 0.0008 0.0084 0.0084 0.0084 可以看到，矩阵A无论是谱半径或是任意范数的值都小于1，可知在 A不变的情况下，Jacobi和Gauss法必然收敛。 2） b取不同的值，x0=20*ones(20,1), e=10^-5, m=50 条 件下的情况对比 迭代 B= B= B= B=20*ones(20,1) B=2000*ones(20,1) 次数 [1:20]’ [10:10:200]’ [20:-1:1]’ Jacobi 24 26 24 23 30 Gauss 16 17 16 15 20 根据1）分析的结果，可以证明无论b取任何值，采用两种方法迭代 均收敛，但b的值的变化会影响迭代的次数；且Gauss迭代法总是比 Jacobi迭代法收敛速度更快。 下表列出的是B=[1:20]’情况下部分结算结果，可以很明显的看到两 种迭代法的收敛速度不同: Jacobi K=1 K=2 K=3 K=4 K=5 K=6 … K=22 K=23 X1 X2 X3 5.3333 9.0000 2.7500 1.5394 1.0874 0.8858 0.7982 … 0.7247 0.7247 4.3333 2.6644 1.9248 1.6082 1.4652 … 1.3444 1.3444 K=2 K=3 K=4 K=5 K=6 … K=14 K=15
由于精度问题，在迭代的最后几次中从显示的数位已经不能看出 标准值与计算值得差别，但是若采用long显示设定，就可以看到更多 位小数的显示，其结果符合最初设定的精度e，数据繁琐，略。 3） x0取不同的值，b=20*ones(20,1), e=10^-5, m=50 条 件下的情况对比 迭代 X0= X0= X0= X0=20*ones(20,1) X0=2000*ones(20,1) 次数 [1:20]’ [10:10:200]’ [20:-1:1]’ Jacobi 22 27 22 23 31 Gauss 15 18 15 15 20 根据1）分析的结果，可以证明无论X0取任何值，采用两种方法迭代 均收敛，但x0的值的变化会影响迭代的次数；且Gauss迭代法总是比 Jacobi迭代法收敛速度更快。 下表列出的是x0=[1:20]’情况下部分结算结果，可以很明显的看到两 种迭代法的收敛速 度不同: Jacobi K=1 K=2 K=3 K=4 K=5 K=6 … K=20 K=21 X1 X2 X3 X1 X2 X3 7.2500 8.6250 9.2228 9.4536 9.5541 9.5974 … 9.6327 9.6327 7.6667 9.9583 10.814 11.183 11.341 11.411 … 11.4683 11.4683 8.1667 10.757 11.816 12.282 12.487 12.579 … 12.6560 12.6560 K=2 K=3 K=4 K=5 K=6 … K=13 K=14 9.6327 7.2500 8.9352 9.4267 9.5720 9.6150 9.6276 … 9.6327
break end end e=eig(bgs) n1=norm(bgs,1); n2=norm(bgs); nn=norm(bgs,inf); min([n1 n2 nn]) 3） 操作函数： %构造矩阵A n=20; a1=sparse(1:n,1:n,3,n,n); 造，比较方便 a2=sparse(1:n-1,2:n,-0.5,n,n); a3=sparse(1:n-2,3:n,-0.25,n,n); a=a1+a2+a3+a2'+a3'; a=full(a);
%按稀疏矩阵的输入法构
%还原为满矩阵
%通过给定不同的初始向量x0或者右端项b，以及规定不同的误差要 求，进行jacobi和gauss %迭代，得到的结果y1、y2位两种迭代的次数，同时输出迭代结果，便 于分析 b= x0= e= m= y1=jacobi(a,b,x0,e,m); y2=gauss(a,b,x0,e,m); 4） 改变矩阵A： 先对jacobi函数作一定修正，方便分析，命名为jacobi2,如下: function y=jacobi2(a,b,x0,e,m) d=diag(diag(a)); u=-triu(a,1); l=-tril(a,-1);
bj=d^-1*(l+u); fj=d^-1*b; x=[x0,zeros(20,m-1)]; n1=norm(bj,1); %计算范数 n2=norm(bj); nn=norm(bj,inf); q=min([n1 n2 nn]); y(1)=q; %输出结果1：范数的最 小值，判断收敛速度的方法 for k=1:m x(:,k+1)=bj*x(:,k)+fj; if norm(x(:,k+1)-x(:,k),inf)<e y(2)=k; %输出结果2：迭 代次数 break end end %改变A矩阵主对角元素的值，比较jacobi迭代的收敛速度，即迭代误 差满足%时的迭代次数 b=(1:20)'; %取定b x0=20*ones(20,1); %取定x0 e=10^-5; %取定精度 r=20; %设置主对角元素 扩大的最大倍数 y=0; m=50; for k=1:r; a1=k*sparse(1:n,1:n,3,n,n); %只需更改A的主对角元 素 a=a1+a2+a3+a2'+a3'; %重新构造A a=full(a); y(k,1:3)=[k,jacobi2(a,b,x0,e,m)]; end y %输出对角线扩大倍 数\最小范数\迭代次数
大倍数 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 9.0000 10.0000 11.0000 12.0000 13.0000 14.0000 15.0000 16.0000 17.0000 18.0000 19.0000 20.0000
%判断迭代后是否满足迭代中止条件： y=x(:,1:k+1); 和迭代结果 sizej=k ; 输出迭代次数 break end 迭代 end
%赋给y所有中间值 %若去掉;号，则 %并结束迭代 %若不成立，继续
%以下部分为验证迭代公式收敛的方法，仅需运行一次即可，因为收敛 性完全由A矩阵决定,而A %在本题是固定不变的；通过判断中B的谱半径或范数大小（B在jacobi 迭代法中%为矩阵bj），即可得知收敛性: e=eig(bj) %输出全部特征值，若，则收敛 n1=norm(bj,1); %计算1-范数 n2=norm(bj); %计算2-范数 nn=norm(bj,inf)； %计算-范数 q=min([n1 n2 nn]) %由于谱半径不超过人以一种范数，所以只要范数的最小值q<1,也可判 断迭代法收敛 2） Gauss迭代法：与Jacobi程序结构相同，不再注释 function y=gauss(a,b,x0,e,m) y=0; d=diag(diag(a)); u=-triu(a,1); l=-tril(a,-1); x=[x0,zeros(20,m-1)]; bgs=(d-l)^-1*u; fgs=(d-l)^-1*b; for k=1:m x(:,k+1)=bgs*x(:,k)+fgs; if norm(x(:,k+1)-x(:,k),inf)<e y=x(:,1:k+1); sizeg=k;
实验 5：线性代数方程组的数值解法
习题3：
已知方程组,其中，定义为： 试通过迭代法求解此方程组，认识迭代法收敛的含义以及迭代初值和方 程组系数矩阵性质对收敛速度的影响。实验要求： （1） 选取不同的初始向量x0和不同的方程组右端向量b，给定迭 代误差要求，用雅可比迭代法和高斯-赛德尔迭代法计算， 观测得到的迭代向量序列是否均收敛？若收敛，记录迭代 次数，分析计算结果并得出结论； （2） 取定右端向量b和初始向量x0，将A的主对角线元素成倍的 增长若干次，非主对角元素不变，每次用雅可比迭代法计 算，要求迭代误差满足,比较收敛速度，分析现象并得出结 论。 1、 程序设计（可直接粘贴运行） 1） Jacobi迭代法 function y=jacobi(a,b,x0,e,m) %定义jacobi函数，其中：a,b为线性方程组中的矩阵和右端向量；x0 为初始值； %e和m分别为人为设定的精度和预计迭代次数；运行结果y为迭代的结 果和所有中间值组成的 %矩阵 y=0; %对y初始化 d=diag(diag(a)); %按雅可比迭代标准形 形式取主对角元素作为矩阵D u=-triu(a,1); %取上三角矩阵u l=-tril(a,-1); %取下三角矩阵l bj=d^-1*(l+u); fj=d^-1*b; x=[x0,zeros(20,m-1)]; %初始化x,其中x1=x0,即 初始值 for k=1:m %人为规定迭代次 数，防止不收敛迭代导致死循环 x(:,k+1)=bj*x(:,k)+fj; %jacobi迭代 if norm(x(:,k+1)-x(:,k),inf)<e