清华大学高等数值计算(李津)实践题目一(共轭梯度CG法,Lanczos算法与MINRES算法)
共轭梯度法步骤
共轭梯度法步骤共轭梯度法是一种求解线性方程组的迭代算法，它以高效稳定的特点而广受欢迎。
以下是共轭梯度法的步骤：步骤1：初始化首先，我们需要有一个初始向量x0和一个初始残量r0=b-Ax0。
其中，A为系数矩阵，b为常数向量。
步骤2：计算方向向量令d0=r0，表示第一次迭代的方向向量。
步骤3：计算步进长度令α0=(r0·r0)/(d0·Ad0)，其中·表示向量的点积。
α0表示迭代过程中每个方向向量的步进长度。
步骤4：更新解向量令x1=x0+α0d0，表示迭代后的解向量。
步骤5：计算新残量令r1=r0-α0Ad0。
步骤6：判断终止条件如果r1的范数小于预设阈值，或者迭代次数达到预设次数，终止迭代。
否则，进入下一次迭代。
步骤7：更新方向向量令β1=(r1·r1)/(r0·r0)，表示更新方向向量的轴线。
步骤8：计算新方向向量令d1=r1+β1d0，表示新的迭代方向向量。
步骤9：计算新的步进长度令α1=(r1·r1)/(d1·Ad1)。
步骤10：更新解向量令x2=x1+α1d1。
步骤11：更新残量令r2=r1-α1Ad1。
步骤12：重复步骤6至11，直至满足终止条件。
总结起来，共轭梯度法的步骤主要包括初始化、计算方向向量、计算步进长度、更新解向量、计算新残量、判断终止条件、更新方向向量、计算新的步进长度、更新解向量和更新残量等。
该算法迭代次数较少，收敛速度快，适用于大规模线性方程组的求解。
共轭梯度算法的设计与实现毕业设计
本科毕业设计（论文）开题报告题目：共轭梯度算法的设计与实现学生姓名：院（系）：理学院专业班级：指导教师：完成时刻：2010年月日共轭梯度算法的设计与实现摘要：共轭梯度法是最优化方式中一种重要的方式.它是介于最速下降法与牛顿法之间的一个方式，它仅需利用一阶导数信息，克服了最速下降法收敛慢的缺点，又幸免了牛顿法需要存储和计算Hesse矩阵并求逆的缺点，共轭梯度法不仅是解决大型线性方程组最有效的方式之一，也是解大型非线性最优化问题最有效的算法之一.最近几年来，随着运算机的飞速进展和实际问题中大规模优化问题的涌现，寻觅快速有效的共轭梯度法成了学者们研究的热点方向之一.本文要紧研究求解无约束最优化问题的共轭梯度法.具体从以下方面着手，介绍所选课题的研究背景、意义和研究现状；论述共轭梯度法理论，包括共轭梯度法的一些大体概念，计算公式和迭代步骤；引入一些实际问题，成立共轭梯度算法；依照已有学者的研究结果，分析算法的收敛性；运用MATLAB编程，代入实例中的相关数据取得一些相关的数值结果；最后对数据进行分析验证，判定该算法的有效性.关键词：无约束最优化；共轭梯度法；迭代；全局收敛性Conjugate Gradient Algorithm Design and ImplementationAbstract：Conjugate gradient method is one of the important optimization methods. It is between the steepest descent method and Newton method, which only requires the first order derivative. It not only overcomes the shortcomings of slow convergence of steepest descent method, but also avoids the shortcomings of storing and calculating the inverse Hesse matrix of Newton method. The conjugate gradient method is not only one of the most useful methods to solve the large-scale linear equations, but also is one of the most effective algorithms to solve a large nonlinear optimization problem.In recent years, with the rapid development of computer and the practical issues in the emergence of large-scale optimization problems, looking for quick and efficient conjugate gradient method becomes one of the most popular directions of scholars. This paper mainly discusses the conjugate gradient method in the unconstrained optimization methods. We focus on the following aspects, describe the background and the significance of the selected research projects; elaborate the theory of conjugate gradient method, which includes some of the basic concepts, formulas and iterative steps; introduce some practical problems, establish the conjugate gradient algorithm; based on the existing academic research results, analyze the convergence of the algorithm; use the MATLAB programming, behalf of the relevant data into the instances to get some relevant numerical results; finally, analyze the data and determine the effectiveness of the algorithm.Key words：Unconstrained optimization；Conjugate gradient method；Iteration；Global convergence目录第一章绪论 ........................................................................................ 错误!未定义书签。
清华大学 李津老师 数值分析第二次实验作业
就不再赘述了。 二、实际计算 生成十个不同的(最好属不同类型或有不同性质的)的 m n 矩阵，这里 m, n 100 ， 用你选择的算法对其做 SVD,比较不同方法的效果(比如计算小气一直和对应左右奇异向量的 误差，效率等)，计算时间和所需存储量等，根据结果提出对算法的认识。 1.误差 在实验中，我们取 m=200，n=100，利用 orth()函数生成了正交矩阵������、������，再生成 了不同奇异值分布的奇异值矩阵������，再通过������ = ������������������，计算出不同的待分解矩阵。 各矩阵奇异值分不如下表所示 序号 1 2 3 4 5 6 7 8 9 10 奇异值个数 10 20 30 40 50 60 70 80 90 100 奇异值分布 10 → 1 20 → 1 30 → 1 40 → 1 50 → 1 60 → 1 70 → 1 80 → 1 90 → 1 100 → 1
−1 −1 −1 −1 −1 −1 −1 −1 −1 −1
经过 matlab 计算，我们得到了两种算法对奇异值的估计误差表，如下所示
序号 svd 1 2 3 4 5 1.2135e-29 1.0336e-28 3.9976e-28 7.9768e-28 1.5711e-27
i 1
100
i
i
2.5232e-27 4.6720e-27 5.8535e-27 8.7958e-27 9.8885e-27
r
i 1
ui uis r
2 2
i 1
vi vis r
2 2
lansvd
svd
lansvd
2 1.6 2.1333 2 2.16 1.8 2.3429 2.15 2 2
清华大学高等数值计算(李津)实践题目二(SVD计算及图像压缩)(包含matlab代码)
主要思想：
传统的 算法和分而治之方法都是基于先将矩阵二对角化的 方法，而 方法是完全不同的方法，它通过一系列平面旋转（也称为 变换）；最终使得矩阵收敛于一个对角矩阵。
其中每一个 设计成，使得对于一对选定的指标 ，有如下式子成立：
最终收敛时，令 ，有 ，且 具有互相正交的列向量，可以将 写成 ，其中 为对角矩阵， 为正交矩阵；从而得到原矩阵的 分解式。
先将矩阵 二对角化得到二对角阵 ，然后将求 的 问题分为两个子问题。
先将矩阵 写成如下形式：
其中 ， 为上二对角矩阵， 是相应维数的向量中的第 个单位向量，并且一般我们取 。
现在假如我们已经求得了 ， 的 如下：
，
并且令 为 的最后一行， 为 的第一行。那么将这些带入 ，可以得到：
可以看出中间的矩阵形式上很简单，除了对角元素和第 行上的元素，其它元素都为 ；这里先把它的 的计算问题放到最后，而假定已知其 为： ，将它带入上式，就可以得到 的 分解式：
主要步骤：
令 是最初的矩阵， 是 经过一次 旋转过得到的矩阵，其中 为对角矩阵， 是列向量范数全为1的矩阵。可以知道对 进行单边的 变换相当于对矩阵 进行双边的 变换。最终得到的矩阵 的列向量的范数就是所要求的奇异值，并且最终归一化以后的矩阵 的列向量就是左奇异向量。
第2部分实验计算
针对以上三种算法用 进行数值模拟，三个算法源自 工具箱，分别为 （传统的SVD算法）， （分而治之）和 （ 方法）。实验中找10幅彩色图片（尺寸为 ），并将其灰度化，转化为数据矩阵，计算3种算法的计算时间和最小奇异值误差，以及左右奇异向量误差，并得到了相应的对比曲线图。左右奇异向量的误差通过 和 的矩阵2-范数的来评估。最小特征值误差，使用默认 自带的SVD算法得到的结果作为参考标准（实际中大多数情况也是无法知道确切的值）。
清华大学高等数值分析作业李津1——矩阵基础
20130917题目求证：在矩阵的LU 分解中，111n n Tn ij i j j i j L I e e α-==+⎛⎫=- ⎪⎝⎭∑∑证明：在高斯消去过程中，假设0jj a ≠ ，若a=0，可以通过列变换使得前面的条件成立，这里不考虑这种情况。
对矩阵A 进行LU 分解，()()()()()1111111L M n M M M n ---=-=••-………… ，其中()1n Tn ij i j i j M j I e e α=+⎛⎫=+ ⎪⎝⎭∑ ，i e 、j e 为n 维线性空间的自然基。
()M j 是通过对单位阵进行初等变换得到，通过逆向的变换则可以得到单位阵，由此很容易得到()M j 的逆矩阵为1n T n ij i j i j I e e α=+⎛⎫- ⎪⎝⎭∑。
故111n n T n ij i j n j i j L I e e I α-==+⎛⎫⎛⎫=- ⎪ ⎪ ⎪⎝⎭⎝⎭∏∑上式中的每一项均是初等变换，从右向左乘，则每乘一次相当于对右边的矩阵进行一次向下乘法叠加的初等变换。
由于最初的矩阵为单位阵，变换从右向左展开，因而每一次变换不改变已经更新的数据，既该变换是从右向左一列一列更新数据，故11nn Tn ij i j j i j L I e e α==+⎛⎫=- ⎪⎝⎭∑∑。
数学证明：1n Tij i j i j e e α=+⎛⎫ ⎪⎝⎭∑具有,000n j j A -⎛⎫ ⎪⎝⎭ 和1,1000n j n j B -+-+⎛⎫⎪⎝⎭ 的形式，且有+1,-11,10000=000n j j n j n j A B --+-+⎛⎫⎛⎫⎪⎪⎝⎭⎝⎭而11n n T ij i j j k i j e e α-==+⎛⎫ ⎪⎝⎭∑∑具有1,1000n k n k B -+-+⎛⎫⎪⎝⎭的形式，因此： 1311111211121==n n n n n n T T T n ij i j n ij i j n ik i k j i j j i j k n i k n n T n i i n ik i i i k L I e e I e e I e e I e e I e ααααα---==+==+=-=+==+⎡⎤⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫=---⎢⎥ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪⎝⎭⎝⎭⎢⎥⎝⎭⎝⎭⎝⎭⎝⎭⎣⎦⎛⎫⎛⎫⎛⎫=-- ⎪ ⎪ ⎝⎭⎝⎝⎭∏∑∏∑∑∑∑∑……11211n n n T Tk n ik i kk k i k e I e e α--===+⎛⎫⎛⎫=- ⎪⎪ ⎪⎭⎝⎭⎝⎭∑∑∑#20130924题目一问：能否用逐次householder 相似变换变实矩阵A 为上三角矩阵，为什么？解：不能用逐次householder 相似变换变A 为上三角矩阵，原因如下：A 记作：()12=,,n A a a a ……， ，存在householder 阵1H s.t. 1111H a e α= ，则()()()111111111111111111111,,,0T Th H AH H a A H e H A H e H A H h H A H ααα⎛⎫'''=== ⎪⎪'⎝⎭⎛⎫''=+ ⎪ ⎪⎝⎭11H A H ''第一列的元素不能保证为1e 的倍数，故无法通过householder 变换实现上三角化。
3_共轭梯度法
共轭梯度法
基本思想
CG的基本思想是把共轭性与最速下降法相结合，利用已 知点处的梯度构造一组共轭方向，并沿着此组方向进行搜索， 求出目标函数的极小点。
d k f ( xk ) k 1d k 1 , k 1,2,..., n 1.
定义--共共轭梯度法的形式 (A) 正定二次函数的无约束最优化问题的共轭梯度法形式
共轭方向法和共轭梯度法
问题1:
如何建立有效的算法？
从二次模型到一般模型.
问题2: 什么样的算法有效呢？ 二次终止性.
简介
共轭方向法和共轭梯度法
共轭梯度法是介于最速下降法与牛顿法之间的一个方法，它仅需利用一阶导 数信息，但克服了最速下降法收敛慢的缺点，又避免了牛顿法需要存储 和计算Hesse矩阵并求逆的缺点，共轭梯度法不仅是解决大型线性方程组 最有用的方法之一，也是解大型非线性最优化最有效的算法之一.
注意
由于Rn中共扼方向最多有 n 个，因此在用上述二种方法求 解目标函数为非二次函数的无约束最优化问题(7.1.1)时，在 n 步之后构造的搜索方向不再是共轭的，从而降低了收敛速度克 服的办法是重设初始点，即把经过 n 次迭代得到的 xn 作为初始 点重新迭代．
算法步骤—FR共轭梯度法
共轭梯度法
Step1: Step2: Step3: Step4: Step5: Step6: Step7:
注： 若 G I ,则
因此, 共轭是正交的推广． 定义--共轭方向法
是正交的，
共轭方向法
性质
特例
n
共轭方向法
算法步骤
Step1: Step2: Step3: Step4:
Step5:
共轭方向法
特例
注
共轭梯度(CG)算法
共轭梯度（CG）算法共轭梯度（Conjugate Gradient, CG）算法是一种用于求解线性方程组的迭代算法。
它主要用于求解对称正定矩阵的线性方程组，如最小二乘问题、PDE（偏微分方程）问题等。
CG算法通过利用矩阵的对称性和正定性，以及向量的共轭关系，实现了高效的求解线性方程组的能力。
CG算法的基本思想是通过一系列共轭的方向，逐步逼近方程组的解。
它利用了矩阵的特性，减少了计算量和存储需求，并且具有较快的收敛速度。
下面将介绍CG算法的原理和过程。
首先，假设我们要求解一个线性方程组Ax=b，其中A是对称正定矩阵，b是已知向量，x是待求解向量。
我们通过迭代的方式逼近x的解，即x(k)。
CG算法的迭代过程如下：1.初始化：选择一个初始解x(0)，设置r(0)=b-Ax(0)，p(0)=r(0)，k=0；2. 迭代计算：计算步长alpha(k)和更新向量x(k+1)：alpha(k) = (r(k)^T * r(k)) / (p(k)^T * A * p(k))x(k+1) = x(k) + alpha(k) * p(k)3. 计算残差向量r(k+1)和比例系数beta(k+1)：r(k+1) = r(k) - alpha(k) * A * p(k)beta(k+1) = (r(k+1)^T * r(k+1)) / (r(k)^T * r(k))4.更新方向p(k+1)：p(k+1) = r(k+1) + beta(k+1) * p(k)5.终止条件判断：如果满足终止条件，停止迭代；否则，令k=k+1，返回步骤2在CG算法中，为了降低数值误差和迭代次数，通常会使用预条件技术，如Jacobi预条件、不完全Cholesky预条件等。
预条件技术可以通过对矩阵进行适当的近似，加速算法的收敛。
CG算法的收敛性和效率主要与矩阵的条件数有关。
对于条件数较大的矩阵，CG算法的迭代次数会增加，收敛速度会减慢。
因此，在实际应用中，通常会选择合适的预条件技术和求解策略，以提高CG算法的效率和稳定性。
共轭梯度实验报告
共轭梯度实验报告共轭梯度实验报告引言：共轭梯度是一种常用的优化算法，广泛应用于数值计算和机器学习等领域。
本实验旨在探究共轭梯度算法的原理和应用，并通过实验验证其在解决线性方程组和最小二乘问题中的有效性和优越性。
一、共轭梯度算法的原理共轭梯度算法是一种迭代法，用于求解对称正定矩阵的线性方程组。
其基本思想是通过选择一组互相共轭的搜索方向，以最小化目标函数的二次型形式。
共轭梯度算法的核心步骤包括初始化、计算搜索方向、计算步长和更新解向量等。
二、共轭梯度算法在线性方程组求解中的应用共轭梯度算法在求解线性方程组方面具有独特的优势。
相比于传统的直接求解方法，共轭梯度算法不需要存储整个矩阵，仅需存储向量和少量中间变量，节省了内存空间。
同时，共轭梯度算法具有较快的收敛速度，能够在有限的迭代次数内得到较精确的解。
三、共轭梯度算法在最小二乘问题中的应用最小二乘问题是一类常见的优化问题，广泛应用于数据拟合和参数估计等领域。
共轭梯度算法在最小二乘问题中的应用主要体现在正规方程法和QR分解法的改进上。
通过共轭梯度算法，可以有效地求解最小二乘问题，得到更准确的拟合结果。
四、实验设计与结果分析本实验选择了一组线性方程组和最小二乘问题进行测试，分别使用共轭梯度算法和传统直接求解方法进行比较。
实验结果表明，共轭梯度算法在求解线性方程组和最小二乘问题时，具有更快的收敛速度和更高的精度。
尤其在大规模问题上，共轭梯度算法的优势更加明显。
结论：共轭梯度算法是一种有效的优化算法，适用于求解对称正定矩阵的线性方程组和最小二乘问题。
通过选择互相共轭的搜索方向，共轭梯度算法能够在有限的迭代次数内得到较精确的解。
在实际应用中，共轭梯度算法具有较快的收敛速度和较高的精度，是一种值得推广和应用的算法。
总结：通过本次实验，我们深入了解了共轭梯度算法的原理和应用，并通过实验验证了其在线性方程组和最小二乘问题中的有效性和优越性。
共轭梯度算法作为一种常用的优化算法，在数值计算和机器学习等领域具有广泛的应用前景。
共轭梯度算法分析与实现
共轭梯度算法分析与实现
梯度下降是一种常用的优化算法，用于求解优化问题。
它通过迭代的
方式不断沿着梯度的反方向更新参数，以最小化损失函数。
然而，梯度下
降算法在处理大规模数据时会变得非常慢，因为它需要计算全部训练样本
的梯度。
为了解决这个问题，共轭梯度算法被提出。
共轭梯度算法是一种适用于解决对称正定矩阵形式下的线性方程组的
优化算法。
它在每一步更新参数时，会按照预先选择好的方向进行更新。
这些方向通常是互相共轭的，这意味着每一个方向都是相对于其他方向来
说是正交的。
共轭梯度算法的原理是，通过每次迭代选择共轭方向来加速
梯度下降算法的收敛速度。
具体而言，共轭梯度算法中的每一步迭代可以分为四个部分：初始化、步长、更新参数和计算残差。
首先，在初始化阶段设定初始参数和初始残差，并选择一个适当的共轭方向。
然后，在步长阶段，通过线方法选择一
个合适的步长。
接下来，在更新参数阶段，根据步长和共轭方向更新参数。
最后，在计算残差阶段，计算新的残差，并检查是否达到停止条件。
如果
没有达到停止条件，那么就继续迭代进行和更新。
共轭梯度算法相对于梯度下降算法有几个优点。
首先，它不需要计算
全部训练样本的梯度，这样可以加速算法的收敛速度。
其次，它可以解决
对称正定矩阵形式下的线性方程组，这在很多实际问题中非常常见。
最后，共轭梯度算法在存储以及计算量上都比较少，所以可以处理大规模数据。
第三章 线性代数方程组的共轭梯度法
对于 k = 0 ,1, 2 ,
tk
=
(r(k) )T p(k) (p(k) )T Ap(k)
x (k+1) = x (k) + t k p(k)
r(k+1) = b − Ax(k+1) = r(k) − t k Ap(k)
sk
=
− (p(k) )T Ar(k+1) (p(k) )T Ap(k)
p(k+1) = r (k+1) + sk p(k)
2：在射线x = x(k ) + tp (k )上求f ( x)的极小点x(k +1) , 即 f ( x(k +1) ) = min f ( x(k ) + tp (k ) )
t>0
3： 判 断 ∇ f ( x (k +1) ) ≤ ε , 或 x(k +1 ) − x(k) ≤ ε ！
是 ！ x* ≈ x (k +1)， 停 止 ; 否 ！ k = k + 1, 转 1.
=
x(0)
+t0p(0).
华北电力大学数理学院
School of mathematics & physics
§3 共轭梯度法
对于 k = 1,2,
现有x ( k )及共轭方向p ( k −1)，则
r(k) = b − Ax(k)
在r ( k )和p ( k −1) 确定的超平面上 找共轭方向
其中
p(k) = r(k) + sk p(k−1)
华北电力大学数理学院
School of mathematics & physics
定理1：设p 1 , p 2 , p n 是关于n阶对称正定 矩阵A共轭的向量组，则以 p(k)为下降方向的算
清华大学高等数值分析(李津)所有作业答案合集
20130917题目求证：在矩阵的LU 分解中，111n n Tn ij i j j i j L I e e α-==+⎛⎫=- ⎪⎝⎭∑∑证明：在高斯消去过程中，假设0jj a ≠ ，若a=0，可以通过列变换使得前面的条件成立，这里不考虑这种情况。
对矩阵A 进行LU 分解，()()()()()1111111L M n M M M n ---=-=∙∙-………… ，其中()1n Tn ij i j i j M j I e e α=+⎛⎫=+ ⎪⎝⎭∑ ，i e 、j e 为n 维线性空间的自然基。
()M j 是通过对单位阵进行初等变换得到，通过逆向的变换则可以得到单位阵，由此很容易得到()M j 的逆矩阵为1n Tn ij i j i j I e e α=+⎛⎫- ⎪⎝⎭∑。
故111n n T n ij i j n j i j L I e e I α-==+⎛⎫⎛⎫=- ⎪ ⎪ ⎪⎝⎭⎝⎭∏∑上式中的每一项均是初等变换，从右向左乘，则每乘一次相当于对右边的矩阵进行一次向下乘法叠加的初等变换。
由于最初的矩阵为单位阵，变换从右向左展开，因而每一次变换不改变已经更新的数据，既该变换是从右向左一列一列更新数据，故11nn Tn ij i j j i j L I e e α==+⎛⎫=- ⎪⎝⎭∑∑。
数学证明：1nTi j i ji j ee α=+⎛⎫ ⎪⎝⎭∑具有,000n j jA -⎛⎫ ⎪⎝⎭ 和1,1000n j n j B -+-+⎛⎫⎪⎝⎭ 的形式，且有+1,-11,10000=000n j j n j n j AB --+-+⎛⎫⎛⎫⎪⎪⎝⎭⎝⎭ 而11n n T ij i j j k i j e e α-==+⎛⎫ ⎪⎝⎭∑∑具有1,1000n k n k B -+-+⎛⎫⎪⎝⎭的形式，因此：1311111211121==n n n n n n T T T n ij i j n ij i j n ik i k j i j j i j k n i k n n T n i i n ik i i i k L I e e I e e I e e I e e I e ααααα---==+==+=-=+==+⎡⎤⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫=---⎢⎥ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪⎝⎭⎝⎭⎢⎥⎝⎭⎝⎭⎝⎭⎝⎭⎣⎦⎛⎫⎛⎫⎛⎫=-- ⎪ ⎪ ⎝⎭⎝⎝⎭∏∑∏∑∑∑∑∑……11211n n n T Tk n ik i kk k i k e I e e α--===+⎛⎫⎛⎫=- ⎪⎪ ⎪⎭⎝⎭⎝⎭∑∑∑#20130924题目一问：能否用逐次householder 相似变换变实矩阵A 为上三角矩阵，为什么？解：不能用逐次householder 相似变换变A 为上三角矩阵，原因如下：A 记作：()12=,,n A a a a ……， ，存在householder 阵1H s.t. 1111H a e α= ，则()()()111111111111111111111,,,0T Th H AH H a A H e H A H e H A H h H A H ααα⎛⎫'''=== ⎪⎪'⎝⎭⎛⎫''=+ ⎪ ⎪⎝⎭11H A H ''第一列的元素不能保证为1e 的倍数，故无法通过householder 变换实现上三角化。
共轭梯度法
一．介绍共轭梯度法（Conjugate Gradient ）是介于最速下降法与牛顿法之间的一个方法，它仅需利用一阶导数信息，但克服了最速下降法收敛慢的缺点，又避免了牛顿法需要存储和计算Hesse 矩阵并求逆的缺点，共轭梯度法不仅是解决大型线性方程组最有用的方法之一，也是解大型非线性最优化最有效的算法之一。
在各种优化算法中，共轭梯度法是非常重要的一种。
其优点是所需存储量小，具有步收敛性，稳定性高，而且不需要任何外来参数。
共轭梯度法中最关键的两点是，搜索方向)(k d 和最佳步长k α。
其基本步骤是在点)(k X 处选取搜索方向)(k d , 使其与前一次的搜索方向)1(-k d 关于A 共轭，即(1)()(1),0k k k d d Ad --<>=然后从点)(k X 出发，沿方向)(k d 求得)(X f 的极小值点)1(+k X , 即)(min )()()(0)1(k dX f X f k k αλ+=>+如此下去, 得到序列{)(k X }。
不难求得0,)1()(>=<-k k Ad d 的解为)()()1(k k k k d X X α+=+其中，><><-=)()()()(,,k k k k kAd d Ad r α注意到)(k d 的选取不唯一，我们可取)1(1)()()(--+-∇=k k k k d X f d β由共轭的定义0,)1()(>=<-k k Ad d 可得：><><-=----)1()1()1()(1,,k k k k k Ad d Ad r β 共轭梯度法的计算公式如下：⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧+=><><-=+=><><-=-=-==+------(k)(k)1)(k )()()()()1(1(k))()1()1()1()(1(k)(k)(0)(0)d X X,,r ,,X r Xr d k k k k k k k k k k k k k k Ad d Ad r d d Ad d Ad r A b A b ααββ 二．程序框图定义矩阵A 和向量bAx=b定义x 的初值将x 代入计算公式误差到达精度要求Yes输出xNo 迭代出新的x 结束开始三．源码n=100;%矩阵阶数，可以按照题目需要更改syms x1 r1 d1A=zeros(n,n);b=zeros(n,1);b(1,1)=-1;b(n,1)=-1;for i=1:nA(i,i)=-2;endfor i=1:n-1A(i,i+1)=1;A(i+1,i)=1;endx1=zeros(n,1);for i=1:n*1000r1=b-A*x1;d1=r1;a=(r1'*d1)/(d1'*A*d1);x1=x1+a*d1;r2=b-A*x1;if(norm(x1)<=eps)breakendbb=-(r2'*A*d1)/(d1'*A*d1);d1=r2+bb*d1;enddisp([x1])四．结果矩阵A100阶的结果200阶的结果400阶的结果。
共轭方向法_实验报告
一、实验目的1. 理解共轭方向法的原理和步骤；2. 掌握共轭方向法在无约束优化问题中的应用；3. 比较共轭方向法与其他优化方法的优缺点。
二、实验原理共轭方向法是一种求解无约束最优化问题的算法，其核心思想是利用共轭方向来迭代搜索最优解。
共轭方向法的主要步骤如下：1. 初始化：选取初始点x0，计算目标函数f(x0)的梯度∇f(x0)；2. 确定搜索方向：计算搜索方向p0=-∇f(x0)；3. 沿搜索方向进行一维搜索：找到一维搜索的最优步长λ0，使得f(x0+λ0p0)≤f(x0)；4. 更新迭代点：计算新的迭代点x1=x0+λ0p0；5. 重复步骤2-4，直到满足终止条件。
共轭方向法要求搜索方向满足共轭条件，即对于给定的对称正定矩阵A，如果两个向量p和q满足(p)Ap=0，则称p和q关于A是共轭的。
三、实验步骤1. 初始化：设置目标函数f(x)，选择初始点x0，计算梯度∇f(x0)；2. 确定搜索方向：计算搜索方向p0=-∇f(x0)；3. 进行一维搜索：使用黄金分割法或二分法找到最优步长λ0；4. 更新迭代点：计算新的迭代点x1=x0+λ0p0；5. 判断终止条件：如果满足|∇f(x1)|<ε或迭代次数超过最大迭代次数，则停止迭代；6. 输出结果：输出最优解x1和对应的目标函数值f(x1)。
四、实验结果与分析1. 实验数据设置目标函数f(x)为f(x)=(x-1)²+(x-2)²，要求找到函数的最小值。
2. 实验结果通过共轭方向法求解，得到最优解x1=1.5，对应的目标函数值f(x1)=0.25。
3. 结果分析共轭方向法在求解该问题时能够找到全局最优解，且迭代次数较少。
与其他优化方法相比，共轭方向法具有以下优缺点：优点：（1）原理简单，易于实现；（2）在多数情况下能够找到全局最优解；（3）迭代次数较少。
缺点：（1）对初始点的选择较为敏感；（2）在迭代过程中可能会出现维度退化现象；（3）当目标函数的梯度变化较大时，搜索效率会降低。
共轭梯度（ConjugateGradient，CG）算法
共轭梯度（ConjugateGradient，CG）算法
以下皆为从⽹络资料获取的感性认知
共轭定义
共轭在数学、物理、化学、地理等学科中都有出现。
本意：两头⽜背上的架⼦称为轭，轭使两头⽜同步⾏⾛。
共轭即为按⼀定的规律相配的⼀对。
通俗点说就是孪⽣。
在数学中有共轭复数、共轭根式、共轭双曲线、共轭矩阵等。
求解⽬标
共轭梯度法是数值分析⾥⾯⼀类⾮常重要的⽅法，是⽤来解决⽆约束凸⼆次规划问题的⼀种⽅法
⽅法的⽬的就是通过迭代求得多元⼆次函数 [公式] 的极值点
基本思想
每个⽅向上只找⼀次，对于 [公式] 维空间就进⾏ [公式] 次计算，也就是说，我每个⽅向上都能⾛到极致，这样我就不会⾛任何的回头路了，效率也就最⾼。
那么是如何实现的呢？如果每个⽅向上都是正交的，我肯定就是在这个⽅向上⾛到极致了，这样⽅向也就确定了，也就是每⼀次都⾛正交⽅向，然后再确定每次的步长就可以了。
抽象图。
lanczos算法及C++实现（一）框架及简单实现
lanczos 算法及C++实现（⼀）框架及简单实现1. lanczos ⽅法的⼤致思路为了求m 阶⽅阵X 最⼤的r 个特征值和特征向量： X m ×m ≈U m ×r S r ×r U T m ×r ,其中U 是列正交矩阵，即 U T U =I ，每⼀列为⼀个特征向量，S 是对⾓阵，对⾓线上每个元素为特征值。
r 为分解的秩lanczos 算法分三步求解:1） 对X 进⾏正交变换得到⼀个三对⾓阵T ：X =PTP T ，其中P 为正交阵T =α0β10000…β1α1β2000…0β2α3β300…00β3α4β40…000β4α5β5…………………… 2） 对T 进⾏奇异值分解：T =W m ×r SW T m ×r3） 最后 X ≈PWS (PW )T 即为所求2. lanczos ⽅法的优点1） T 矩阵是⼀个三对⾓阵，很稀疏，因此对它的特征值分解会⾮常快。
2） 如果r 很⼩，可以不⽤求整个T ，⽽是求其左上1.5r 阶的⽅阵即可得到很好的近似。
这样对T 的特征值分解会更快。
() 另外⼀种⽅法动态决定求多少lanczos 向量， 这⾥简单总结⼀下，初始选取n ≤m 个lanczos 向量，然后求T 的最⼤特征值和最⼩特征值，然后逐渐增加n ，并更新这两个特征值，直到这种更新⾮常⼩为⽌。
如果⽤在求r 个主特征向量，那么可以更新第1个和第r 个特征值。
3. lanczos ⽅法的导出考虑幂迭代产⽣的⼀系列向量P =[v ,Xv ,X 2v ,…X m −1v ]，其中v 是任意向量。
假设X 的特征值分解为X =∑i λi ξi ξT i ，并且特征值按绝对值从⼤到⼩排序 |λ1|≥|λ2|≥…。
则 Xv =∑i λi ξi ξT i v =∑i λi v ′i ξi ，其中v ′i =ξT i v 是v 在正交基[ξ1,ξ2…]下的坐标。
共轭梯度求解器的FPGA设计与实现
共轭梯度求解器的FPGA设计与实现宋庆增;顾军华【期刊名称】《计算机应用》【年(卷),期】2011(031)009【摘要】To overcome the disadvantage of inefficient and bad real-time capability in software version Conjugate Gradient (CG) iterative solver, a CC iterative solver was designed and implemented on Field Programmable Gate Array ( FPGA) platform. The design of CG iterative solver was based on hardware/software co-design. Hardware CG co-processor implemented the code of enormous computation and simple control, which could accelerate the system. The code of control complexity and less calculation was still performed in the microprocessor. The use of row-interleaved data flow could make the system not to stall to improve performance. The experimental results illustrate that hardware CG iterative solver can speed up about 5. 7 times over the software version of the same algorithm.%针时共轭梯度(CG)迭代算法软件执行效率低、实时性差的缺点,提出一种基于现场可编程逻辑门阵列(FPGA)平台的CG迭代求解器.设计采用软硬件结合的方式构建整个系统,CG协处理器执行CG迭代算法中计算量大、控制简单的代码,以达到硬件加速的目的.控制复杂、计算量较少的代码则依旧在微处理上执行.设计采用行交错数据流,使得整个系统完全无停顿的运行,提高了计算性能.实验结果表明,与软件执行相比,硬件CG协处理器可以获得最高5.7倍的性能加速.【总页数】4页(P2571-2573,2588)【作者】宋庆增;顾军华【作者单位】河北工业大学电气工程学院,天津300401;河北工业大学计算机科学与软件学院,天津300401【正文语种】中文【中图分类】TP31【相关文献】1.基于Web Services的开放式分布计算求解器的设计与实现 [J], 何伟;孔梦荣;车战斌2.基于可满足性问题求解器的星上FPGA永久损伤容错技术研究 [J], 孙兆伟;刘源;赵丹;陈健;张世杰3.基于FPGA的Jacobi迭代求解器研究 [J], 宋庆增;顾军华;张金珠4.基于不完全算法的并行FPGA SAT求解器 [J], 黎铁军;马柯帆;张建民5.基于不完全分解预优共轭梯度法的电源和地线网络求解器 [J], 武晓海;殷莉;洪先龙因版权原因，仅展示原文概要，查看原文内容请购买。
共轭梯度法解油藏数值模拟线性方程中的稀疏矩阵及C++程序(测试通过)
for(i=0;i<=n-1;i++) {
x[i]=0.0; p[i]=b[i]; r[i]=b[i]; } i=0; while(i<=n-1) { for(k=0;k<=n-1;k++) {
s[k]=0.0; for(j=0;j<=n-1;j++)
s[k]=s[k]+a[k*n+j]*p[j]; } d=0.0; e=0.0; for(k=0;k<=n-1;k++) {
Байду номын сангаасfor(int jj=0;jj<5;jj++) {
aa[m++]=a[ii][jj];
} } double b[5]={23,32,33,31,66}; eps=0.000001; grad(5,aa,b,eps,x); for(i=0;i<=4;i++) {
cout<<"x["<<i<<"]="<<x[i]<<endl; } }
d=d+p[k]*b[k]; e=e+p[k]*s[k]; } alpha=d/e; for(k=0;k<=n-1;k++) { x[k]=x[k]+alpha*p[k];
} for(k=0;k<=n-1;k++) {
q[k]=0.0; for(j=0;j<=n-1;j++) {
q[k]=q[k]+a[k*n+j]*x[j]; } } d=0.0; for(k=0;k<=n-1;k++) { r[k]=b[k]-q[k]; d=d+r[k]*s[k]; } beta=d/e; d=0.0; for(k=0;k<=n-1;k++) { d=d+r[k]*r[k]; } d=sqrt(d); if(d<eps) {
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
高等数值计算实践题目一1. 实践目的本次计算实践主要是在掌握共轭梯度法，Lanczos 算法与MINRES 算法的基础上，进一步探讨这3种算法的数值性质，主要研究特征值特征向量对算法收敛性的影响。
2. 实践过程（一）生成矩阵（1）作5个100阶对角阵i D 如下：1D 对角元：1,1,...,20,1+0.1(-20),21,...,100j j d j d j j ====2D 对角元：1,1,...,20,1+(-20),21,...,100j j d j d j j ==== 3D 对角元：,1,...,80,81,81,...,100j j d j j d j ====4D 对角元：,1,...,40,41,41,...,60,41+(60),61,...,100j j j d j j d j d j j =====-= 5D 对角元：,1,...,100j d j j ==记i D 的最大模特征值和最小模特征值分别为1i λ和i n λ，则i D 特征值分布有如下特点：1D 的特征值有较多接近于i n λ，并且1/i i n λλ较小，2D 的特征值有较多接近于i n λ，并且1/i i n λλ较大， 3D 的特征值有较多接近于1i λ，并且1/i i n λλ较大，4D 的特征值有较多接近于中间模特征值，并且1/i i n λλ较大， 5D 的特征值均匀分布，并且1/i i n λλ较大（2）随机生成10个100阶矩阵j M ：(100(100))j M fix rand =并作它们的QR 分解，得j Q 和j R ，这样可得50个对称的矩阵Tij j i j A Q DQ =，其中i D 的对角元就是ij A 的特征值，若它们都大于0，则ij A 正定，j Q 的列就是相应的特征向量。
结合（1）可知，ij A 都是对称正定阵。
（二）计算结果以下计算，均选定精确解(100,1)exact x ones =，初值0(100,1)x zeros =由ij exact kA x b =计算得到k b （算法中要求解的精度为10e -）。
Lanczos 算法和MINRES 算法均分别使用不带m 循环和有带m 循环进行计算，其中带m 循环时，取15m =进行计算。
（1） 特征值分布对3种算法的影响特征值分布对共轭梯度法的影响特征值分布对()Lanczos m 算法的影响特征值分布对()MINRES m 算法的影响（2） 特征向量矩阵对3种算法的影响特征向量对共轭梯度法的影响（取定2D ）特征向量对Lanczos 算法的影响（取定2D ）特征向量对()Lanczos m 算法的影响（取定2D ）特征向量对MINRES 算法的影响（取定2D ）特征向量对()MINRES m 算法的影响（取定2D ）（3）从ij A 选取5个画出收敛率曲线 1.共轭梯度法的1k A k Ae k e -曲线以及kAe k 曲线（当迭代步数()j j k n n >表示在步收敛，1k A kAk Ae e e -和由于初始化均为零）nczos 算法的1k A k Ae k e -曲线以及kAe k 曲线（当迭代步数()j j k n n >表示在步收敛，1k A kAk Ae e e -和由于初始化均为零）3.()Lanczos m 算法的1k A k Ak e -曲线以及kAe k 曲线（当迭代步数()j j k n n >表示在步收敛，1k A kAk Ae e e -和由于初始化均为零）4.MINRES 算法的1k A k Ak e -曲线以及kAe k 曲线（当迭代步数()j j k n n >表示在步收敛，1k A kAk Ae e e -和由于初始化均为零）5.()MINRES m 算法的1k A k Ae k e -曲线以及kAe k 曲线（当迭代步数()j j k n n >表示在步收敛，1k A kAk Ae e e -和由于初始化均为零）（三）相关讨论（1）由矩阵i D 的生成方法，可以知道：1D 的特征值有较多接近于i n λ，并且1/i i n λλ较小，2D 的特征值有较多接近于i n λ，并且1/i i n λλ较大，3D 的特征值有较多接近于1i λ，并且1/i i n λλ较大，4D 的特征值有较多接近于中间模特征值，并且1/i i n λλ较大，5D 的特征值均匀分布，并且1/i i n λλ较大。
观察分析计算结果，有以下几点结论：1.当1/i i n λλ较小时，3种算法都收敛的特别快；当1/i i n λλ较大时，3种算法的收敛速度都要慢一些，这与理论分析是一致的。
2.当特征值有较多接近i n λ与特征值有较多接近于1i λ以及特征值有较多接近中间值时，几乎不影响3种算法的收敛速度，而且特征向量也几乎不影响3种算法的收敛速度。
（2）1.在填特征值分布对3种算法的影响的表中，对于ij A 的选择原则：1.1首先由于要研究特征值分布对3种算法的影响，因此选取的ij A 的特征值应该照顾的不同的特征值分布；1.2对于固定的特征值分布i D 对应的ij A 的生成方法可知，由于j M 是随机的，j Q 也是随机的，因此，选取了1j =的情况，即1121314151,,,,A A A A A 。
2.在填特征值向量对3种算法的影响的表中，取定一个i D 的选择原则：考虑矩阵i D 的特征值分布以及对计算结果的观察（特征值向量对3种算法的影响不大），于是取定2D 。
3.在画收敛率曲线图时，对于ij A 的选择原则：同1的选取原则（3）当矩阵阶数n （本次计算100n =）不是很大时，对于对称正定阵，3种算法的收敛速度都比较快，尤其是当1/i i n λλ较小时，3种算法都收敛的特别快，这与课程上的理论分析是一致的。
其次，若选取合适的m （比如本次计算所选取的=15m ），使用Lanczos m ()算法和MINRES m ()算法都有相当令人满意的收敛性质，而且相比共轭梯度法的收敛速度要快得多。
另外，也可以看到共轭梯度法数值上是十分不稳定的，相比较而言，使用Lanczos m ()算法和MINRES m ()算法要稳定的多。
（4）stationary 迭代法中的SOR 以及使用对称的SOR 公式构造的预处理CG 法的迭代结果如下所示：（算法中要求解的精度为5e -）观察分析计算结果：stationary 迭代法中的SOR 的收敛速度要慢很多，而且解的精度不是很高，同时数值上也不是很稳定；使用预处理CG 法，可以加快收敛速度，并且可以增加数值稳定性。
因此，相比较而言，使用预处理CG 法与Lanczos m ()算法和MINRES m ()算法要好得多。
stationary 迭代法中的SOR 算法SOR CG1.stationary 迭代法中的SOR 算法（当迭代步数()j j k n n >表示在步收敛，1k A kAk Ae e e -和由于初始化均为零）2.使用对称的SOR 公式构造的预处理CG 法（当迭代步数()j j k n n >表示在步收敛，1k A kAk Ae e e -和由于初始化均为零）（5）对于非正定的对称正定系数矩阵，比较Lanczos m ()算法和MINRES m ()算法 1.构造非正定的对称系数矩阵，类似正定对称矩阵的构造方式，选取对角阵6D ，对角元：1,1,...,20,(1+(-20)),21,...,100j j d j d j j ===-=，随机矩阵11(100(100))M fix rand =，作它的QR 分解，得11Q ，令11611T AA Q D Q =，则AA 为非正定的对称阵。
2.选定精确解(100,1)exact x ones =，初值0(100,1)x zeros =由ij exact k A x b =计算得到k b （算法中要求解的精度为10e -）。
L a n c z o s 算法和MINRES 算法均分别使用不带m 循环和有带m 循环进行计算，其中带m 循环时，取15m =进行计算。
3.由于系数矩阵非正定，,Axx Ax =<>会出现负值，故不能作为范数，此时考虑使用2-范数来代替，画出Lanczos 算法和MINRES 算法的212k k e k e -曲线以及2ke k 曲线。
4.结合图表可以看出，对于此非正定的对称系数矩阵AA ，2种算法均有很好的收敛性质，不仅收敛速度特别快，而且解的精度可以很高。
同时，可以看出，Lanczos 算法的收敛曲线图在开始时出现了一个较大的峰值，相应的误差突然增大，表现出数值不稳定。
因此，对于本例，MINRES 算法要略优于Lanczos 算法，MINRES 算法要相对稳定得多。
（当迭代步数()j j k n n >表示在步收敛，2212k k k e e -和由于初始化均为零）（当迭代步数()j j k n n >表示在步收敛，2212k k k e e -和由于初始化均为零）DMAQ.m 用于生成矩阵D,M,A,QCG.m 共轭梯度法Lanczos1.m Lanczos算法Lanczos.m Lanczos(m)算法MINRES1.m MINRES算法MINRES.m MINRES(m)算法SOR.m SOR算法SOR_CG.m 对称SOR的预处理共轭梯度法main.m 主函数，包括画图代码NPDS.m 用于非正定的对称系数矩阵的算例Lanczos1_AA.m 用于非正定的对称系数矩阵的Lanczos算法Lanczos_AA.m 用于非正定的对称系数矩阵的Lanczos(m)算法MINRES1_AA.m 用于非正定的对称系数矩阵的MINRES算法MINRES_AA.m 用于非正定的对称系数矩阵的MINRES(m)算法function [D,Q,A]=DMAQa12=linspace(1.1,9,80);d1=[a11,a12];D1=diag(d1);a21=ones(1,20);a22=linspace(2,81,80);d2=[a21,a22];D2=diag(d2);a31=linspace(1,80,80);a32=81*ones(1,20);d3=[a31,a32];D3=diag(d3);a41=linspace(1,40,40); a42=41*ones(1,20);a43=linspace(42,81,40); d4=[a41,a42,a43];D4=diag(d4);d5=1:100;D5=diag(d5);N=100;M1=fix(100*rand(N)); [Q1,R1]=qr(M1);M2=fix(100*rand(N)); [Q2,R2]=qr(M2);M3=fix(100*rand(N)); [Q3,R3]=qr(M3);M4=fix(100*rand(N)); [Q4,R4]=qr(M4);M5=fix(100*rand(N)); [Q5,R5]=qr(M5);M6=fix(100*rand(N)); [Q6,R6]=qr(M6);M7=fix(100*rand(N)); [Q7,R7]=qr(M7);M8=fix(100*rand(N)); [Q8,R8]=qr(M8);M9=fix(100*rand(N)); [Q9,R9]=qr(M9);M10=fix(100*rand(N)); [Q10,R10]=qr(M10);D=cat(3,D1,D2,D3,D4,D5);Q=cat(3,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,Q10);k=0;for i=1:5for j=1:10k=k+1;A(:,:,k)=Q(:,:,j)*D(:,:,i)*Q(:,:,j)';endendendfunction [x,iter,w1,w2]=CG(A,b,x0,xx)iter=0;tol=1.0e-10;[N1,N2]=size(A);w1=zeros(N1,1);w2=zeros(N1,1);x=x0;r=b-A*x;p=r;for i=1:N1t1=x-xx;alpha=(r'*r)/(p'*(A*p));r1=r;x=x+alpha*p;r=r-alpha*A*p;belta=(r'*r)/(r1'*r1);p=r+belta*p;t2=x-xx;w1(i,1)=sqrt((t2'*(A*t2))/(t1'*(A*t1))); w2(i,1)=sqrt(t2'*(A*t2));iter=iter+1;if norm(t2)<tolbreak;endendfunction [x,iter,w1,w2]=Lanczos1(A,bb,x0,xx)iter=0;tol=1.0e-10;[N1,N2]=size(A);m=N1;w1=zeros(N1,1);w2=zeros(N1,1);x=x0;q=zeros(N1);a=zeros(1,m);b=zeros(1,m);r0=bb-A*x0;b0=norm(r0);q(:,1)=r0/b0;a(1)=q(:,1)'*A*q(:,1);r1=A*q(:,1)-a(1)*q(:,1);b(1)=norm(r1);q(:,2)=r1/b(1);for i=2:N1a(i)=q(:,i)'*A*q(:,i);r1=A*q(:,i)-a(i)*q(:,i)-b(i-1)*q(:,i-1); b(i)=norm(r1);b(i)=b(i);q(:,i+1)=r1/b(i);m=i;Q=q(:,1:m);H=Q'*A*Q;y=H\(Q'*r0);z=Q*y;t1=x-xx;x=x0+z;t2=x-xx;w1(i,1)=sqrt((t2'*(A*t2))/(t1'*(A*t1))); w2(i,1)=sqrt(t2'*(A*t2));iter=iter+1;if norm(t2)<tolbreak;endendendfunction [x,iter,w1,w2]=Lanczos(A,bb,x0,xx)iter=0;tol=1.0e-10;[N1,N2]=size(A);w1=zeros(N1,1);w2=zeros(N1,1);x=x0;for j=1:N1m=15;q=zeros(N1);a=zeros(1,m);b=zeros(1,m);r0=bb-A*x;b0=norm(r0);q(:,1)=r0/b0;a(1)=q(:,1)'*A*q(:,1);r1=A*q(:,1)-a(1)*q(:,1);b(1)=norm(r1);q(:,2)=r1/b(1);for i=2:ma(i)=q(:,i)'*A*q(:,i);r1=A*q(:,i)-a(i)*q(:,i)-b(i-1)*q(:,i-1); b(i)=norm(r1);b(i)=b(i);q(:,i+1)=r1/b(i);endQ=q(:,1:m);H=Q'*A*Q;y=H\(Q'*r0);z=Q*y;t1=x-xx;x=x0+z;t2=x-xx;w1(j,1)=sqrt((t2'*(A*t2))/(t1'*(A*t1)));w2(j,1)=sqrt(t2'*(A*t2));iter=iter+1;if norm(t2)<tolbreak;endx0=x;endendfunction [x,iter,w1,w2]=MINRES1(A,bb,x0,xx)iter=0;tol=1.0e-10;[N1,N2]=size(A);w1=zeros(N1,1);w2=zeros(N1,1);x=x0;m=N1;q=zeros(N1);a=zeros(1,m);b=zeros(1,m);r0=bb-A*x0;b0=norm(r0);q(:,1)=r0/b0;a(1)=q(:,1)'*A*q(:,1);r1=A*q(:,1)-a(1)*q(:,1);b(1)=norm(r1);q(:,2)=r1/b(1);for i=2:N1a(i)=q(:,i)'*A*q(:,i);r1=A*q(:,i)-a(i)*q(:,i)-b(i-1)*q(:,i-1); b(i)=norm(r1);b(i)=b(i);q(:,i+1)=r1/b(i);m=i;Q=q(:,1:m);HH=zeros(m+1,m);for j=1:mHH(j,j)=a(j);endfor j=1:m-1HH(j,j+1)=b(j);endfor j=2:mHH(j,j-1)=b(j-1);endHH(m+1,m)=b(m);H=HH;[Q1,R1]=qr(H);RR=R1(1:m,:);e1=zeros(m+1,1);e1(1,1)=1;bb1=b0*Q1'*e1;bb2=bb1(1:m,1);y=RR\(bb2);z=Q*y;t1=x-xx;x=x0+z;t2=x-xx;w1(i,1)=sqrt((t2'*(A*t2))/(t1'*(A*t1))); w2(i,1)=sqrt(t2'*(A*t2));iter=iter+1;if norm(t2)<tolbreak;endendendfunction [x,iter,w1,w2]=MINRES(A,bb,x0,xx)iter=0;tol=1.0e-10;[N1,N2]=size(A);w1=zeros(N1,1);w2=zeros(N1,1);x=x0;for j=1:N1m=15;q=zeros(N1);a=zeros(1,m);b=zeros(1,m);r0=bb-A*x;b0=norm(r0);q(:,1)=r0/b0;a(1)=q(:,1)'*A*q(:,1);r1=A*q(:,1)-a(1)*q(:,1);b(1)=norm(r1);q(:,2)=r1/b(1);for i=2:ma(i)=q(:,i)'*A*q(:,i);r1=A*q(:,i)-a(i)*q(:,i)-b(i-1)*q(:,i-1); b(i)=norm(r1);b(i)=b(i);q(:,i+1)=r1/b(i);endQ=q(:,1:m);QQ=q(:,1:(m+1));H=QQ'*A*Q;[Q1,R1]=qr(H);RR=R1(1:m,:);e1=zeros(m+1,1);e1(1,1)=1;bb1=(b0*Q1'*e1)';bb2=bb1(1,1:m);y=RR\(bb2');z=Q*y;t1=x-xx;x=x0+z;t2=x-xx;w1(j,1)=sqrt((t2'*(A*t2))/(t1'*(A*t1)));w2(j,1)=sqrt(t2'*(A*t2));iter=iter+1;if norm(t2)<tolbreak;endx0=x;endendfunction [x,iter,w1,w2]=SOR(A,b,x0,xx)iter=0;tol=1.0e-4;w=1.5;[N1,N2]=size(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-w*L)\((1-w)*D+w*U);f=(D-w*L)\b*w;x=x0;w1=zeros(N1,1);w2=zeros(N1,1);for i=1:N1t1=x-xx;x=B*x+f;t2=x-xx;iter=iter+1;if norm(t2)<tolbreak;endw1(i,1)=sqrt((t2'*(A*t2))/(t1'*(A*t1))); w2(i,1)=sqrt(t2'*(A*t2));endendfunction [x,iter,w1,w2]=SOR_CG(A,b,x0,xx)iter=0;%µü´ú²½Êýtol=1.0e-4;%Îó²î¿ØÖÆ[N1,N2]=size(A);w1=zeros(N1,1);w2=zeros(N1,1);w=1.2;D=diag(diag(A));L=-tril(A,-1);M=(1/(w*(2-w)))*(D-w*L)*(D^(-1))*(D-w*L');x=x0;r=b-A*x;z=M\r;p=z;for i=1:N1alpha=(r'*z)/(p'*A*p);t1=x-xx;x=x+alpha*p;r1=r-alpha*A*p;z1=M\r1;belta=(z1'*z1)/(r'*z);p=z1+belta*p;r=r1;z=z1;t2=x-xx;iter=iter+1;if norm(t2)<tolbreak;endw1(i,1)=sqrt((t2'*(A*t2))/(t1'*(A*t1))); w2(i,1)=sqrt(t2'*(A*t2));endendfunction [x,iter,w1,w2]=SOR_CG(A,b,x0,xx)iter=0;tol=1.0e-4;[N1,N2]=size(A);w1=zeros(N1,1);w2=zeros(N1,1);w=1.2;D=diag(diag(A));L=-tril(A,-1);M=(1/(w*(2-w)))*(D-w*L)*(D^(-1))*(D-w*L');x=x0;r=b-A*x;z=M\r;p=z;for i=1:N1alpha=(r'*z)/(p'*A*p);t1=x-xx;x=x+alpha*p;r1=r-alpha*A*p;z1=M\r1;belta=(z1'*z1)/(r'*z);p=z1+belta*p;r=r1;z=z1;t2=x-xx;iter=iter+1;if norm(t2)<tolbreak;endw1(i,1)=sqrt((t2'*(A*t2))/(t1'*(A*t1))); w2(i,1)=sqrt(t2'*(A*t2));endendclc;clear;[D,Q,A]=DMAQ;xx=ones(100,1);x0=zeros(100,1);B=zeros(100,50);KK=zeros(50,2);X_CG=zeros(100,50);Iter_CG=zeros(50,1);W1_CG=zeros(100,50);W2_CG=zeros(100,50);EE_CG=zeros(50,2);X_Lanczos1=zeros(100,50);Iter_Lanczos1=zeros(50,1);W1_Lanczos1=zeros(100,50);W2_Lanczos1=zeros(100,50);EE_Lanczos1=zeros(50,2);X_Lanczos=zeros(100,50);Iter_Lanczos=zeros(50,1);W1_Lanczos=zeros(100,50);W2_Lanczos=zeros(100,50);EE_Lanczos=zeros(50,2);X_MINRES1=zeros(100,50);Iter_MINRES1=zeros(50,1);W1_MINRES1=zeros(100,50);W2_MINRES1=zeros(100,50);EE_MINRES1=zeros(50,2);X_MINRES=zeros(100,50);Iter_MINRES=zeros(50,1);W1_MINRES=zeros(100,50);W2_MINRES=zeros(100,50);EE_MINRES=zeros(50,2);X_SOR=zeros(100,50);Iter_SOR=zeros(50,1);W1_SOR=zeros(100,50);W2_SOR=zeros(100,50);EE_SOR=zeros(50,2);X_SOR_CG=zeros(100,50);Iter_SOR_CG=zeros(50,1);W1_SOR_CG=zeros(100,50);W2_SOR_CG=zeros(100,50);EE_SOR_CG=zeros(50,2);for k=1:50K=A(:,:,k);temp=cond(K);KK(k,1)=temp;KK(k,2)=(temp-1)/(temp+1);B(:,k)=A(:,:,k)*xx;[X_CG(:,k),Iter_CG(k,1),W1_CG(:,k),W2_CG(:,k)]=CG(K,B(:,k),x0,xx); WW=W1_CG(:,k);EE_CG(k,1)=mean(WW(1:Iter_CG(k,1)));EE_CG(k,2)=norm(X_CG(:,k)-xx);[X_Lanczos(:,k),Iter_Lanczos(k,1),W1_Lanczos(:,k),W2_Lanczos(:,k)]=La nczos(K,B(:,k),x0,xx);WW=W1_Lanczos(:,k);EE_Lanczos(k,1)=mean(WW(1:Iter_Lanczos(k,1)));EE_Lanczos(k,2)=norm(X_Lanczos(:,k)-xx);[X_Lanczos1(:,k),Iter_Lanczos1(k,1),W1_Lanczos1(:,k),W2_Lanczos1(:,k) ]=Lanczos1(K,B(:,k),x0,xx);WW=W1_Lanczos1(:,k);EE_Lanczos1(k,1)=mean(WW(1:Iter_Lanczos1(k,1)));EE_Lanczos1(k,2)=norm(X_Lanczos1(:,k)-xx);[X_MINRES(:,k),Iter_MINRES(k,1),W1_MINRES(:,k),W2_MINRES(:,k)]=MINRES (K,B(:,k),x0,xx);WW=W1_MINRES(:,k);EE_MINRES(k,1)=mean(WW(1:Iter_MINRES(k,1)));EE_MINRES(k,2)=norm(X_MINRES(:,k)-xx);[X_MINRES1(:,k),Iter_MINRES1(k,1),W1_MINRES1(:,k),W2_MINRES1(:,k)]=MI NRES1(K,B(:,k),x0,xx);WW=W1_MINRES1(:,k);EE_MINRES1(k,1)=mean(WW(1:Iter_MINRES1(k,1)));EE_MINRES1(k,2)=norm(X_MINRES1(:,k)-xx);[X_SOR(:,k),Iter_SOR(k,1),W1_SOR(:,k),W2_SOR(:,k)]=SOR(K,B(:,k),x0,xx );WW=W1_SOR(:,k);EE_SOR(k,1)=mean(WW(1:Iter_SOR(k,1)));EE_SOR(k,2)=norm(X_SOR(:,k)-xx);[X_SOR_CG(:,k),Iter_SOR_CG(k,1),W1_SOR_CG(:,k),W2_SOR_CG(:,k)]=SOR_CG(K,B(:,k),x0,xx);WW=W1_SOR_CG(:,k);EE_SOR_CG(k,1)=mean(WW(1:Iter_SOR_CG(k,1)));EE_SOR_CG(k,2)=norm(X_SOR_CG(:,k)-xx);end¨subplot(211);MN=100;TT=1:MN;plot(TT,W1_CG(1:MN,1),TT,W1_CG(1:MN,11),TT,W1_CG(1:MN,21),TT,W1_CG(1: MN,31),TT,W1_CG(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('{||e_{k}||_{A}}/{||e_{k-1}||_{A}}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('¹²éîÌݶȷ¨');subplot(212);MN=100;TT=1:MN;plot(TT,W2_CG(1:MN,1),TT,W2_CG(1:MN,11),TT,W2_CG(1:MN,21),TT,W2_CG(1: MN,31),TT,W2_CG(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('||e_{k}||_{A}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('¹²éîÌݶȷ¨');subplot(211);MN=100;TT=1:MN;plot(TT,W1_Lanczos1(1:MN,1),TT,W1_Lanczos1(1:MN,11),TT,W1_Lanczos1(1: MN,21),TT,W1_Lanczos1(1:MN,31),TT,W1_Lanczos1(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('{||e_{k}||_{A}}/{||e_{k-1}||_{A}}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('LanczosËã·¨');subplot(212);MN=100;TT=1:MN;plot(TT,W2_Lanczos1(1:MN,1),TT,W2_Lanczos1(1:MN,11),TT,W2_Lanczos1(1: MN,21),TT,W2_Lanczos1(1:MN,31),TT,W2_Lanczos1(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('||e_{k}||_{A}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('LanczosËã·¨');subplot(211);MN=100;TT=1:MN;plot(TT,W1_Lanczos(1:MN,1),TT,W1_Lanczos(1:MN,11),TT,W1_Lanczos(1:MN, 21),TT,W1_Lanczos(1:MN,31),TT,W1_Lanczos(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('{||e_{k}||_{A}}/{||e_{k-1}||_{A}}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('Lanczos(m)Ëã·¨');subplot(212);MN=100;TT=1:MN;plot(TT,W2_Lanczos(1:MN,1),TT,W2_Lanczos(1:MN,11),TT,W2_Lanczos(1:MN, 21),TT,W2_Lanczos(1:MN,31),TT,W2_Lanczos(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('||e_{k}||_{A}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('Lanczos(m)Ëã·¨');subplot(211);MN=100;TT=1:MN;plot(TT,W1_MINRES1(1:MN,1),TT,W1_MINRES1(1:MN,11),TT,W1_MINRES1(1:MN, 21),TT,W1_MINRES1(1:MN,31),TT,W1_MINRES1(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('{||e_{k}||_{A}}/{||e_{k-1}||_{A}}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('MINRESËã·¨');subplot(212);MN=100;TT=1:MN;plot(TT,W2_MINRES1(1:MN,1),TT,W2_MINRES1(1:MN,11),TT,W2_MINRES1(1:MN, 21),TT,W2_MINRES1(1:MN,31),TT,W2_MINRES1(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('||e_{k}||_{A}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('MINRESËã·¨');subplot(211);MN=100;TT=1:MN;plot(TT,W1_MINRES(1:MN,1),TT,W1_MINRES(1:MN,11),TT,W1_MINRES(1:MN,21) ,TT,W1_MINRES(1:MN,31),TT,W1_MINRES(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('{||e_{k}||_{A}}/{||e_{k-1}||_{A}}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('MINRES(m)Ëã·¨');subplot(212);MN=100;TT=1:MN;plot(TT,W2_MINRES(1:MN,1),TT,W2_MINRES(1:MN,11),TT,W2_MINRES(1:MN,21) ,TT,W2_MINRES(1:MN,31),TT,W2_MINRES(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('||e_{k}||_{A}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('MINRES(m)Ëã·¨');subplot(211);MN=100;TT=1:MN;plot(TT,W1_SOR(1:MN,1),TT,W1_SOR(1:MN,11),TT,W1_SOR(1:MN,21),TT,W1_SO R(1:MN,31),TT,W1_SOR(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('{||e_{k}||_{A}}/{||e_{k-1}||_{A}}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('SORËã·¨');subplot(212);MN=100;TT=1:MN;plot(TT,W2_SOR(1:MN,1),TT,W2_SOR(1:MN,11),TT,W2_SOR(1:MN,21),TT,W2_SO R(1:MN,31),TT,W2_SOR(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('||e_{k}||_{A}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('SORËã·¨');subplot(211);MN=100;TT=1:MN;plot(TT,W1_SOR_CG(1:MN,1),TT,W1_SOR_CG(1:MN,11),TT,W1_SOR_CG(1:MN,21),TT,W1_SOR_CG(1:MN,31),TT,W1_SOR_CG(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('{||e_{k}||_{A}}/{||e_{k-1}||_{A}}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('¶Ô³ÆSORµÄÔ¤´¦Àí¹²éîÌݶȷ¨');subplot(212);MN=100;TT=1:MN;plot(TT,W2_SOR_CG(1:MN,1),TT,W2_SOR_CG(1:MN,11),TT,W2_SOR_CG(1:MN,21) ,TT,W2_SOR_CG(1:MN,31),TT,W2_SOR_CG(1:MN,41));xlabel('µü´ú²½Êýk');ylabel('||e_{k}||_{A}');legend('A_{11}','A_{21}','A_{31}','A_{41}','A_{51}');title('¶Ô³ÆSORµÄÔ¤´¦Àí¹²éîÌݶȷ¨');clc;clear;N=100;M11=fix(100*rand(N));[Q11,R11]=qr(M11);a61=ones(1,20);a62=-linspace(2,81,80);d6=[a61,a62];D6=diag(d6);AA=Q11*D6*Q11';xx=ones(100,1);x0=zeros(100,1);bb=AA*xx;norm_AA=zeros(2,1);norm1_AA=zeros(2,1);[X_Lanczos_AA,Iter_Lanczos_AA,W1_Lanczos_AA,W2_Lanczos_AA]=Lanczos_AA (AA,bb,x0,xx);norm_AA(1,1)=norm(X_Lanczos_AA-xx);[X_Lanczos1_AA,Iter_Lanczos1_AA,W1_Lanczos1_AA,W2_Lanczos1_AA]=Lanczo s1_AA(AA,bb,x0,xx);norm1_AA(1,1)=norm(X_Lanczos1_AA-xx);[X_MINRES_AA,Iter_MINRES_AA,W1_MINRES_AA,W2_MINRES_AA]=MINRES_AA(AA,b b,x0,xx);norm_AA(2,1)=norm(X_MINRES_AA-xx);[X_MINRES1_AA,Iter_MINRES1_AA,W1_MINRES1_AA,W2_MINRES1_AA]=MINRES1_AA (AA,bb,x0,xx);norm1_AA(2,1)=norm(X_MINRES1_AA-xx);%»-ͼMN=100;TT=1:MN;subplot(211);plot(TT,W1_Lanczos1_AA,TT,W1_MINRES1_AA)xlabel('µü´ú²½Êýk');ylabel('{||e_{k}||_{2}}/{||e_{k-1}||_{2}}');legend('LanczosËã·¨','MINRESËã·¨');title('LanczosËã·¨ÓëMINRESËã·¨');subplot(212);plot(TT,W2_Lanczos1_AA,TT,W2_MINRES1_AA)xlabel('µü´ú²½Êýk');ylabel('||e_{k}||_{2}');legend('LanczosËã·¨','MINRESËã·¨');title('LanczosËã·¨ÓëMINRESËã·¨');subplot(211);plot(TT,W1_Lanczos_AA,TT,W1_MINRES_AA)xlabel('µü´ú²½Êýk');ylabel('{||e_{k}||_{2}}/{||e_{k-1}||_{2}}');legend('Lanczos(m)Ëã·¨','MINRES(m)Ëã·¨');title('Lanczos(m)Ëã·¨ÓëMINRES(m)Ëã·¨');subplot(212);plot(TT,W2_Lanczos_AA,TT,W2_MINRES_AA)xlabel('µü´ú²½Êýk');ylabel('||e_{k}||_{2}');legend('Lanczos(m)Ëã·¨','MINRES(m)Ëã·¨');title('Lanczos(m)Ëã·¨ÓëMINRES(m)Ëã·¨');function [x,iter,w1,w2]=Lanczos1_AA(A,bb,x0,xx)iter=0;tol=1.0e-10;[N1,N2]=size(A);m=N1;w1=zeros(N1,1);w2=zeros(N1,1);x=x0;q=zeros(N1);a=zeros(1,m);b=zeros(1,m);r0=bb-A*x0;b0=norm(r0);q(:,1)=r0/b0;a(1)=q(:,1)'*A*q(:,1);r1=A*q(:,1)-a(1)*q(:,1);b(1)=norm(r1);q(:,2)=r1/b(1);for i=2:N1a(i)=q(:,i)'*A*q(:,i);r1=A*q(:,i)-a(i)*q(:,i)-b(i-1)*q(:,i-1);b(i)=norm(r1);b(i)=b(i);q(:,i+1)=r1/b(i);m=i;Q=q(:,1:m);H=Q'*A*Q;y=H\(Q'*r0);z=Q*y;t1=x-xx;x=x0+z;t2=x-xx;w1(i,1)=norm(t2)/norm(t1);w2(i,1)=norm(t2);iter=iter+1;if norm(t2)<tolbreak;endendendfunction [x,iter,w1,w2]=Lanczos_AA(A,bb,x0,xx)iter=0;tol=1.0e-10;[N1,N2]=size(A);w1=zeros(N1,1);w2=zeros(N1,1);x=x0;for j=1:N1m=15;q=zeros(N1);a=zeros(1,m);b=zeros(1,m);r0=bb-A*x;b0=norm(r0);q(:,1)=r0/b0;a(1)=q(:,1)'*A*q(:,1);r1=A*q(:,1)-a(1)*q(:,1);b(1)=norm(r1);q(:,2)=r1/b(1);for i=2:ma(i)=q(:,i)'*A*q(:,i);r1=A*q(:,i)-a(i)*q(:,i)-b(i-1)*q(:,i-1);b(i)=norm(r1);b(i)=b(i);q(:,i+1)=r1/b(i);endQ=q(:,1:m);H=Q'*A*Q;y=H\(Q'*r0);z=Q*y;t1=x-xx;x=x0+z;t2=x-xx;w1(j,1)=norm(t2)/norm(t1);w2(j,1)=norm(t2);iter=iter+1;if norm(t2)<tolbreak;endx0=x;endendfunction [x,iter,w1,w2]=MINRES1_AA(A,bb,x0,xx)iter=0;tol=1.0e-10;[N1,N2]=size(A);w1=zeros(N1,1);w2=zeros(N1,1);x=x0;m=N1;q=zeros(N1);a=zeros(1,m);b=zeros(1,m);r0=bb-A*x0;b0=norm(r0);q(:,1)=r0/b0;a(1)=q(:,1)'*A*q(:,1);r1=A*q(:,1)-a(1)*q(:,1);b(1)=norm(r1);q(:,2)=r1/b(1);for i=2:N1a(i)=q(:,i)'*A*q(:,i);r1=A*q(:,i)-a(i)*q(:,i)-b(i-1)*q(:,i-1);b(i)=norm(r1);b(i)=b(i);q(:,i+1)=r1/b(i);m=i;Q=q(:,1:m);HH=zeros(m+1,m);for j=1:mHH(j,j)=a(j);endfor j=1:m-1HH(j,j+1)=b(j);endfor j=2:mHH(j,j-1)=b(j-1);endHH(m+1,m)=b(m);H=HH;[Q1,R1]=qr(H);RR=R1(1:m,:);e1=zeros(m+1,1);e1(1,1)=1;bb1=b0*Q1'*e1;bb2=bb1(1:m,1);y=RR\(bb2);z=Q*y;t1=x-xx;x=x0+z;t2=x-xx;w1(i,1)=norm(t2)/norm(t1);w2(i,1)=norm(t2);iter=iter+1;if norm(t2)<tolbreak;endendendfunction [x,iter,w1,w2]=MINRES_AA(A,bb,x0,xx)iter=0;tol=1.0e-10;[N1,N2]=size(A);w1=zeros(N1,1);w2=zeros(N1,1);x=x0;for j=1:N1m=15;q=zeros(N1);a=zeros(1,m);b=zeros(1,m);r0=bb-A*x;b0=norm(r0);q(:,1)=r0/b0;a(1)=q(:,1)'*A*q(:,1);r1=A*q(:,1)-a(1)*q(:,1);b(1)=norm(r1);q(:,2)=r1/b(1);for i=2:ma(i)=q(:,i)'*A*q(:,i);r1=A*q(:,i)-a(i)*q(:,i)-b(i-1)*q(:,i-1); b(i)=norm(r1);b(i)=b(i);q(:,i+1)=r1/b(i);endQ=q(:,1:m);QQ=q(:,1:(m+1));H=QQ'*A*Q;[Q1,R1]=qr(H);RR=R1(1:m,:);e1=zeros(m+1,1);e1(1,1)=1;bb1=(b0*Q1'*e1)';bb2=bb1(1,1:m);y=RR\(bb2');z=Q*y;t1=x-xx;x=x0+z;t2=x-xx;w1(j,1)=norm(t2)/norm(t1);w2(j,1)=norm(t2);iter=iter+1;if norm(t2)<tolbreak;endx0=x;endend。