下载文档原格式
共轭梯度法步骤共轭梯度法是一种求解线性方程组的迭代算法，它以高效稳定的特点而广受欢迎。
以下是共轭梯度法的步骤：步骤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目录第一章绪论 ........................................................................................ 错误!未定义书签。
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 变换实现上三角化。
共轭梯度（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分解法的改进上。
通过共轭梯度算法，可以有效地求解最小二乘问题，得到更准确的拟合结果。
四、实验设计与结果分析本实验选择了一组线性方程组和最小二乘问题进行测试，分别使用共轭梯度算法和传统直接求解方法进行比较。
实验结果表明，共轭梯度算法在求解线性方程组和最小二乘问题时，具有更快的收敛速度和更高的精度。
尤其在大规模问题上，共轭梯度算法的优势更加明显。
结论：共轭梯度算法是一种有效的优化算法，适用于求解对称正定矩阵的线性方程组和最小二乘问题。
通过选择互相共轭的搜索方向，共轭梯度算法能够在有限的迭代次数内得到较精确的解。
在实际应用中，共轭梯度算法具有较快的收敛速度和较高的精度，是一种值得推广和应用的算法。
总结：通过本次实验，我们深入了解了共轭梯度算法的原理和应用，并通过实验验证了其在线性方程组和最小二乘问题中的有效性和优越性。
共轭梯度算法作为一种常用的优化算法，在数值计算和机器学习等领域具有广泛的应用前景。
共轭梯度算法分析与实现
梯度下降是一种常用的优化算法，用于求解优化问题。
它通过迭代的
方式不断沿着梯度的反方向更新参数，以最小化损失函数。
然而，梯度下
降算法在处理大规模数据时会变得非常慢，因为它需要计算全部训练样本
的梯度。
为了解决这个问题，共轭梯度算法被提出。
共轭梯度算法是一种适用于解决对称正定矩阵形式下的线性方程组的
优化算法。
它在每一步更新参数时，会按照预先选择好的方向进行更新。
这些方向通常是互相共轭的，这意味着每一个方向都是相对于其他方向来
说是正交的。
共轭梯度算法的原理是，通过每次迭代选择共轭方向来加速
梯度下降算法的收敛速度。
具体而言，共轭梯度算法中的每一步迭代可以分为四个部分：初始化、步长、更新参数和计算残差。
首先，在初始化阶段设定初始参数和初始残差，并选择一个适当的共轭方向。
然后，在步长阶段，通过线方法选择一
个合适的步长。
接下来，在更新参数阶段，根据步长和共轭方向更新参数。
最后，在计算残差阶段，计算新的残差，并检查是否达到停止条件。
如果
没有达到停止条件，那么就继续迭代进行和更新。
共轭梯度算法相对于梯度下降算法有几个优点。
首先，它不需要计算
全部训练样本的梯度，这样可以加速算法的收敛速度。
其次，它可以解决
对称正定矩阵形式下的线性方程组，这在很多实际问题中非常常见。
最后，共轭梯度算法在存储以及计算量上都比较少，所以可以处理大规模数据。
高等数值计算实践题目一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。
