代数方程组数值解法

线性代数方程组的数值解法讨论解线性方程组的方法，主要分为直接方法和迭代方法两种。

直接法是在没有舍入误差的假设下能在预定的运算次数内求得精确解。

而实际上，原始数据的误差和运算的舍入误差是不可以避免的，实际上获得的也是近似解。

迭代法是构造一定的递推格式，产生逼近精确解的序列。

对于高阶方程组，如一些偏微分方程数值求解中出现的方程组，采用直接法计算代价比较高，迭代法则简单又实用，因此比较受工程人员青睐。

小组成员本着工程应用，讨论将学习的理论知识转变为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 分解法等，另一类是迭代法（近似法），包括雅克比迭代法、高斯迭代法、超松弛迭代法等。

Hilbert 代数方程组的数值解法1 Hilbert 矩阵的条件数和矩阵的阶数的关系编制Matlab 程序：clear,clc format long %format short for n=1:14; A=hilb(n); %A=rand(n);condA(n)=cond(A,inf); endplot(log10(condA)) 得出图形图1：Hilbert 矩阵和阶次的关系由图可见a.Hilbert 矩阵的2-条件数的对数与阶次近似正比;b.阶次超过12后，计算困难，舍入误差会导致计算不准 结论：Hilbert 矩阵的2-条件数随阶次指数增长，关系大概是： Log 10(Cond(A))= 1.33791720780254(n-1)0246810121424681012141618nl o g 10(c o n d A )条件数的对数-阶次2.各种求解方法的对比2.1 直接法：Gauss 消去方法（程序见附录，下同）在双精度型变量下，即eps=2.220446049250313e-016，对于直接法Gauss 消去方法，求解Hilbert 矩阵方程组，在阶次n=13时，解得: X=[1 0.99997 1.0013 0.97881 1.1906 -0.035441 4.6171 -7.3967 14.089 -12.542 9.9162 -2.3817 1.5624]出现明显错误，可见Gauss 消去方法对病态问题比较敏感。

求解阶次不高。

2.2 Jacobi 迭代法很遗憾，jakobi 迭代法的迭代矩阵的谱半径随阶次线性上升（程序见附录2），普遍大于1，计算结果发送，无法迭代出满意的结果。

需放弃图2：Jacobi 迭代矩阵的谱半径2.3 Gauss-Seidel 迭代与SOR 迭代求解先分析SOR 迭代矩阵的谱半径和松弛因子w 的关系0510152025510152025阶次谱半径Jacobi 迭代收敛性分析图3：SOR 迭代矩阵的谱半径和松弛因子w 的关系可以发现SOR 迭代和G-S 迭代的谱半径都普遍接近于1，因此松弛因子的选择对计算结果的影响不大。

线性代数方程组的数值解法【实验目的】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消去法在消去过程中选择模最大的主元能够得到比较精确的解。

数学中的代数数值方法代数数值方法是数学中一类重要的计算方法，它主要用于解决代数方程、代数方程组和代数函数的数值计算问题。

代数数值方法融合了代数与数值计算的思想，既能通过代数技巧对问题进行化简，又能利用数值计算的方法获得近似解。

一、代数数值方法的概述代数数值方法主要包括以下几个方面：方程求根、代数方程组的数值解法、插值与逼近以及数值微积分等。

在实际问题的计算中，这些方法常用于求解无法或难以用解析方法求解的方程和函数的数值解，以及近似求解函数的值、函数的导数等。

二、方程求根的代数数值方法方程求根是代数数值方法的一个重要分支，常用的方法有二分法、牛顿法和割线法等。

1. 二分法：二分法是一种简单且有效的求解方程根的方法。

它通过迭代逼近的方式不断缩小方程根所在的区间，直至达到预设的精度要求。

2. 牛顿法：牛顿法基于泰勒级数展开和迭代逼近的思想，通过对原函数进行逼近的线性插值，从而求得函数根的近似解。

3. 割线法：割线法与牛顿法类似，但其使用的迭代形式为割线逼近，更加灵活，并且对初值选取的要求相对较低。

三、代数方程组的数值解法代数方程组的数值解法主要包括直接法和迭代法两种。

1. 直接法：直接法是通过消元和替代等操作，将原方程组转化为等价的简化形式，从而获得方程组的解。

常见的直接法包括高斯消元法、克拉默法则等。

2. 迭代法：迭代法通过设置初始迭代值，利用逐步逼近的思想求解方程组。

常见的迭代法有雅可比迭代法、高斯-赛德尔迭代法等。

四、插值与逼近插值与逼近是代数数值方法中的重要内容，它们可以用于生成一些函数的近似模型，从而方便进一步的数值计算。

1. 插值：插值方法可以通过已知数据点之间的插值多项式来逼近实际函数。

常见的插值方法有拉格朗日插值、牛顿插值等。

2. 逼近：逼近方法则是通过已知的离散数据点，寻求逼近函数与原函数之间的最佳拟合。

最小二乘逼近是一种常用的逼近方法。

五、数值微积分数值微积分是代数数值方法的另一个重要分支，它主要用于计算函数的导数、积分以及微分方程的数值解。

线性代数方程组数值解法及MATLAB 实现综述廖淑芳 20122090 数计学院 12计算机科学与技术1班（职教本科） 一、分析课题随着科学技术的发展，提出了大量复杂的数值计算问题，在建立电子计算机成为数值计算的主要工具以后，它以数字计算机求解数学问题的理论和方法为研究对象。

其数值计算中线性代数方程的求解问题就广泛应用于各种工程技术方面。

因此在各种数据处理中，线性代数方程组的求解是最常见的问题之一。

关于线性代数方程组的数值解法一般分为两大类：直接法和迭代法。

直接法就是经过有限步算术运算，可求的线性方程组精确解的方法（若计算过程没有舍入误差），但实际犹如舍入误差的存在和影响，这种方法也只能求得近似解，这类方法是解低阶稠密矩阵方程组级某些大型稀疏矩阵方程组的有效方法。

直接法包括高斯消元法，矩阵三角分解法、追赶法、平方根法。

迭代法就是利用某种极限过程去逐步逼近线性方程组精确解的方法。

迭代法具有需要计算机的存储单元少，程序设计简单，原始系数矩阵在计算过程始终不变等优点，但存在收敛性级收敛速度问题。

迭代法是解大型稀疏矩阵方程组（尤其是微分方程离散后得到的大型方程组）的重要方法。

迭代法包括Jacobi 法SOR 法、SSOR 法等多种方法。

二、研究课题-线性代数方程组数值解法 一、 直接法 1、 Gauss 消元法通过一系列的加减消元运算，也就是代数中的加减消去法，以使A 对角线以下的元素化为零，将方程组化为上三角矩阵；然后，再逐一回代求解出x 向量。

1.1消元过程1. 高斯消元法(加减消元)：首先将A 化为上三角阵，再回代求解。

11121121222212n n n n nn n a a a b a a a b a a a b ⎛⎫ ⎪ ⎪ ⎪ ⎪⎝⎭L L M M O M M L (1)(1)(1)(1)(1)11121311(2)(2)(2)(2)222322(3)(3)(3)3333()()000000nn n n n nn n a a a a b a a a b a a b a b ⎛⎫ ⎪ ⎪ ⎪ ⎪⎪ ⎪⎝⎭L L L M M M OM M L 步骤如下：第一步：1111,2,,i a i i n a -⨯+=L 第行第行 11121121222212n nn n nn n a a a b a a a b a a a b ⎛⎫⎪ ⎪ ⎪⎪⎝⎭L L M M O M M L111211(2)(2)(2)2222(2)(2)(2)200n nn nn n a a a b a a b a a b ⎛⎫⎪⎪ ⎪ ⎪⎝⎭LL M M O M M L第二步：(2)2(2)222,3,,i a i i n a -⨯+=L 第行第行111211(2)(2)(2)2222(2)(2)(2)200nnn nn n a a a b a a b a a b ⎛⎫ ⎪⎪ ⎪ ⎪⎝⎭L L M M O M M L11121311(2)(2)(2)(2)222322(3)(3)(3)3333(3)(3)(3)300000n n nn nn n a a a a b a a a b a a b a a b ⎛⎫⎪ ⎪ ⎪⎪⎪ ⎪⎝⎭LL LM M M O M M L 类似的做下去，我们有：第k 步：()()k ,1,,k ikk kka i i k n a -⨯+=+L 第行第行。

实验5 线性代数方程组的数值解法分1 黄浩 43一、实验目的1.学会用MATLAB软件数值求解线性代数方程组，对迭代法的收敛性和解的稳定性作初步分析；2.通过实例学习用线性代数方程组解决简化的实际问题。

二、For personal use only in study and research; not forcommercial use三、四、实验内容1.《数学实验》第二版（问题1）问题叙述：通过求解线性方程组，理解条件数的意义和方程组性态对解的影响，其中是n阶范德蒙矩阵，即是n阶希尔伯特矩阵，b1，b2分别是的行和。

（1）编程构造（可直接用命令产生）和b1，b2；你能预先知道方程组和的解吗？令n=5，用左除命令求解（用预先知道的解可验证程序）。

（2）令n=5,7,9，…，计算和的条件数。

为观察他们是否病态，做以下试验：b1，b2不变，和的元素，分别加扰动后求解；和不变，b1，b2的分量b1(n),b2(n)分别加扰动后求解。

分析A与b的微小扰动对解的影响。

取10^-10,10^-8,10^-6。

（3）经扰动得到的解记做，计算误差，与用条件数估计的误差相比较。

模型转换及实验过程：（1）小题.由b1,b2为,的行和，可知方程组和的精确解均为n 行全1的列向量。

在n=5的情况下，用matlab编程（程序见四.1），构造,和b1，b2，使用高斯消去法得到的解x1,x2及其相对误差e1,e2（使用excel计算而得）为：由上表可见，当n=5时，所得的解都接近真值，误差在10^-12的量级左右。

（2）小题分别取n=5,7,9,11,13,15，计算和的条件数c1和c2，（程序见四.2），结果如下：由上表可见，二者的条件数都比较大，可能是病态的。

为证实和是否为病态，先保持b不变，对做扰动，得到该情况下的高斯消元解，（程序见四.3），结果如下：(为使结果清晰简洁，在此仅列出n=5,9,13的情况，n=7,11,15略去)=10^-10时：=10^-8时：=10^-6时：由上表可见：a)对于希尔伯特阵，随着阶数的增加，微小扰动对解带来的影响越来越大，到了n=9时，已经有了6倍误差的解，到了n=13时，甚至出现了22倍误差的解元素；而随着的增加，解的偏差似乎也有增加的趋势，但仅凭上述表格无法具体判断（在下一小题中具体叙述）。