MATLAB中cgs函数用法

已于 2024-08-18 10:52:22 修改
阅读量411 收藏 8
点赞数 16
分类专栏： Matlab 文章标签： matlab 算法 开发语言
于 2024-08-18 10:46:11 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/jk_101/article/details/133164491
版权

目录

语法

说明

示例

线性系统的迭代解

使用指定了预条件子的 cgs

提供初始估计值

使用函数句柄代替数值矩阵

        cgs函数的功能是求解线性系统 - 共轭梯度二乘法。

语法

x = cgs(A,b)
x = cgs(A,b,tol)
x = cgs(A,b,tol,maxit)
x = cgs(A,b,tol,maxit,M)
x = cgs(A,b,tol,maxit,M1,M2)
x = cgs(A,b,tol,maxit,M1,M2,x0)
[x,flag] = cgs(___)
[x,flag,relres] = cgs(___)
[x,flag,relres,iter] = cgs(___)
[x,flag,relres,iter,resvec] = cgs(___)

说明

        x = cgs(A,b) 尝试使用共轭梯度二乘法求解关于 x 的线性系统 A*x = b。如果尝试成功，cgs 会显示一条消息来确认收敛。如果 cgs 无法在达到最大迭代次数后收敛或出于任何原因暂停，则会显示一条包含相对残差 norm(b-A*x)/norm(b) 以及该方法停止时的迭代次数的诊断消息。

        x = cgs(A,b,tol) 指定该方法的容差。默认容差是 1e-6。

        x = cgs(A,b,tol,maxit) 指定要使用的最大迭代次数。如果 cgs 无法在 maxit 次迭代内收敛，将显示诊断消息。

        x = cgs(A,b,tol,maxit,M) 指定预条件子矩阵 M 并通过有效求解关于 y 的方程组 AM^−1y=b 来计算 x，其中 y=Mx。使用预条件子矩阵可以改善问题的数值属性和计算的效率。

        x = cgs(A,b,tol,maxit,M1,M2) 指定预条件子矩阵 M 的因子，使得 M = M1*M2。

        x = cgs(A,b,tol,maxit,M1,M2,x0) 指定解向量 x 的初始估计值。默认值为由零组成的向量。

        [x,flag] = cgs(___) 返回一个标志，指示算法是否成功收敛。当 flag = 0 时，收敛成功。您可以将此输出语法用于之前的任何输入参数组合。如果指定了 flag 输出，cgs 将不会显示任何诊断消息。

        [x,flag,relres] = cgs(___) 还会返回相对残差 norm(b-A*x)/norm(b)。如果 flag 为 0，则 relres <= tol。

        [x,flag,relres,iter] = cgs(___) 还会返回计算出 x 时的迭代次数 iter。

        [x,flag,relres,iter,resvec] = cgs(___) 还会在每次迭代中返回残差范数向量（包括第一个残差 norm(b-A*x0)）。

示例

线性系统的迭代解

        使用采用默认设置的 cgs 求解系数矩阵为方阵的线性系统，然后在求解过程中调整使用的容差和迭代次数。

        创建密度为 50% 的随机稀疏矩阵 A。另为 Ax=b 的右侧创建随机向量。

rng default
A = sprand(600,600,.5);
A = A'*A + speye(size(A));
b = rand(600,1);

        使用 cgs 求解 Ax=b。输出显示包括相对残差 的值。

x = cgs(A,b);
cgs stopped at iteration 20 without converging to the desired tolerance 1e-06
because the maximum number of iterations was reached.
The iterate returned (number 20) has relative residual 0.068.

        默认情况下，cgs 使用 20 次迭代和容差 1e-6，对于此矩阵，算法无法在 20 次迭代后收敛。由于残差仍然很大，这说明需要更多的迭代（或预条件子矩阵）。您也可以使用更大的容差，使算法更容易收敛。

使用容差 1e-4 和 40 次迭代再次求解方程组。

x = cgs(A,b,1e-4,40);
cgs stopped at iteration 40 without converging to the desired tolerance 0.0001
because the maximum number of iterations was reached.
The iterate returned (number 40) has relative residual 0.0024.

        即使有更宽松的容差和更多迭代，cgs 也不会收敛。当迭代算法以这种方式停滞时，显然需要预条件子矩阵。

计算 A 的不完全 Cholesky 分解，并使用 L 因子作为 cgs 的预条件子输入。

L = ichol(A);
x = cgs(A,b,1e-4,40,L);
cgs converged at iteration 14 to a solution with relative residual 1.4e-05.

        使用预条件子可以充分改进问题的数值属性，使 cgs 能够收敛。

使用指定了预条件子的 cgs

        检查使用指定了预条件子矩阵的 cgs 来求解线性系统的效果。加载 west0479，它是一个非对称的 479×479 实稀疏矩阵。

load west0479
A = west0479;

定义 b 以使 Ax=b 的实际解是全为 1 的向量。

b = sum(A,2);

设置容差和最大迭代次数。

tol = 1e-12;
maxit = 20;

使用 cgs 根据请求的容差和迭代次数求解。指定五个输出以返回有关求解过程的信息：

  • x 是计算 A*x = b 所得的解。

  • fl0 是指示算法是否收敛的标志。

  • rr0 是计算的解 x 的相对残差。

  • it0 是计算 x 时所用的迭代次数。

  • rv0 是 ‖b−Ax‖ 的残差历史记录组成的向量。

[x,fl0,rr0,it0,rv0] = cgs(A,b,tol,maxit);
fl0
fl0 = 1
rr0
rr0 = 1
it0
it0 = 0

        cgs 未在请求的 20 次迭代内收敛至请求的容差 1e-12，因此 fl0 为 1。实际上，cgs 的行为太差，因此初始估计值 x0 = zeros(size(A,2),1) 是最佳解，并会返回最佳解（如 it0 = 0 所示）。

        为了有助于缓慢收敛，可以指定预条件子矩阵。由于 A 是非对称的，请使用 ilu 生成预条件子 M=L U。指定调降容差，以忽略值小于 1e-6 的非对角线元。通过指定 L 和 U 作为 cgs 的输入，求解预条件方程组 A M^−1M x=b。

setup = struct('type','ilutp','droptol',1e-6);
[L,U] = ilu(A,setup);
[x1,fl1,rr1,it1,rv1] = cgs(A,b,tol,maxit,L,U);
fl1
fl1 = 0
rr1
rr1 = 4.3850e-14
it1
it1 = 3

        在第三次迭代中，使用 ilu 预条件子产生的相对残差小于规定的容差 1e-12。输出 rv1(1) 为 norm(b)，输出 rv1(end) 为 norm(b-A*x1)。

        可以通过绘制每次迭代的相对残差来跟踪 cgs 的进度。绘制每个解的残差历史记录图，并添加一条表示指定容差的线。

semilogy(0:length(rv0)-1,rv0/norm(b),'-o')
hold on
semilogy(0:length(rv1)-1,rv1/norm(b),'-o')
yline(tol,'r--');
legend('No preconditioner','ILU preconditioner','Tolerance','Location','East')
xlabel('Iteration number')
ylabel('Relative residual')

如图所示：

提供初始估计值

        检查向 cgs 提供解的初始估计值的效果。

        创建一个三对角稀疏矩阵。使用每行的总和作为 Ax=b 右侧的向量，使 x 的预期解是由 1 组成的向量。

n = 900;
e = ones(n,1);
A = spdiags([e 2*e e],-1:1,n,n);
b = sum(A,2);

        使用 cgs 求解 Ax=b 两次：一次是使用默认的初始估计值，一次是使用解的良好初始估计值。对两次求解均使用 200 次迭代和默认容差。将第二种求解中的初始估计值指定为所有元素都等于 0.99 的向量。

maxit = 200;
x1 = cgs(A,b,[],maxit);
cgs converged at iteration 17 to a solution with relative residual 8.8e-07.
x0 = 0.99*e;
x2 = cgs(A,b,[],maxit,[],[],x0);
cgs converged at iteration 4 to a solution with relative residual 8e-07.

        在这种情况下，提供初始估计值可以使 cgs 更快地收敛。

返回中间结果

        还可以通过在 for 循环中调用 cgs 来使用初始估计值获得中间结果。每次调用求解器都会执行几次迭代，并存储计算出的解。然后，将该解用作下一批迭代的初始向量。

        例如，以下代码会循环执行四次，每次执行 100 次迭代，并在 for 循环中每通过一次后均存储解向量：

x0 = zeros(size(A,2),1);
tol = 1e-8;
maxit = 100;
for k = 1:4
    [x,flag,relres] = cgs(A,b,tol,maxit,[],[],x0);
    X(:,k) = x;
    R(k) = relres;
    x0 = x;
end

        X(:,k) 是在 for 循环的第 k 次迭代时计算的解向量，R(k) 是该解的相对残差。

使用函数句柄代替数值矩阵

        通过为 cgs 提供用来计算 A*x 的函数句柄（而非系数矩阵 A）来求解线性系统。

        gallery 生成的 Wilkinson 测试矩阵之一是 21×21 三对角矩阵。预览该矩阵。

A = gallery('wilk',21)
A = 21×21

    10     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     1     9     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     1     8     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     1     7     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     1     6     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     1     5     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     1     4     1     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     1     3     1     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     1     2     1     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     1     1     1     0     0     0     0     0     0     0     0     0     0
      ⋮

        Wilkinson 矩阵有特殊的结构，因此您可以用函数句柄来表示 A*x 运算。当 A 乘以向量时，所得向量中的大多数元素为零。结果中的非零元素对应于 A 的非零三对角元素。此外，只有主对角线具有不等于 1 的非零值。

表达式 Ax 变为：

如图所示：

结果向量可以写为三个向量的和：

        在 MATLAB® 中，编写一个函数来创建这些向量并将它们相加，从而给出 A*x 的值：

function y = afun(x)
y = [0; x(1:20)] + ...
    [(10:-1:0)'; (1:10)'].*x + ...
    [x(2:21); 0];
end

（该函数作为局部函数保存在示例的末尾。）

        现在，通过为 cgs 提供用于计算 A*x 的函数句柄，求解线性系统 Ax=b。使用容差 1e-12 和 50 次迭代。

b = ones(21,1);
tol = 1e-12;  
maxit = 50;
x1 = cgs(@afun,b,tol,maxit)
cgs converged at iteration 11 to a solution with relative residual 1.3e-14.
x1 = 21×1

    0.0910
    0.0899
    0.0999
    0.1109
    0.1241
    0.1443
    0.1544
    0.2383
    0.1309
    0.5000
      ⋮

        检查 afun(x1) 是否产生由 1 组成的向量。

afun(x1)
ans = 21×1

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
      ⋮

局部函数

function y = afun(x)
y = [0; x(1:20)] + ...
    [(10:-1:0)'; (1:10)'].*x + ...
    [x(2:21); 0];
end

共轭梯度二乘法

        开发共轭梯度二乘 (CGS) 算法是为了改进双共轭梯度 (BiCG) 算法。CGS 算法不使用残差及其共轭，而是通过使用残差平方来避免使用系数矩阵的转置 [1]。​

        在计算成本相当的情况下，CGS 的收敛速度快于 BiCG，但可能具有不规则的收敛行为，尤其是当初始估计值接近解时 [1]。

​提示

  • ​大多数迭代方法的收敛取决于系数矩阵的条件数 cond(A)。当 A 是方阵时，您可以使用 equilibrate 来改进其条件数，它本身就能使大多数迭代求解器更容易收敛。但如果随后会对经平衡处理的矩阵 B = R*P*A*C 进行因式分解，使用 equilibrate 还可以获得质量更好的预条件子矩阵。

  • 可以使用矩阵重新排序函数（如 dissect 和 symrcm）来置换系数矩阵的行和列，并在系数矩阵被分解以生成预条件子时最小化非零值的数量。这可以减少后续求解预条件线性系统所需的内存和时间。

参考

        [1] Barrett, R., M. Berry, T. F. Chan, et al., Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, SIAM, Philadelphia, 1994.

        [2] Sonneveld, Peter, “CGS: A fast Lanczos-type solver for nonsymmetric linear systems,” SIAM J. Sci. Stat. Comput., January 1989, Vol. 10, No. 1, pp. 36–52.

jk_101
关注
  • 16
    点赞
  • 踩
  • 8
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
matlab中cgs收敛曲线,求解线性方程组 - 共轭梯度二乘法
weixin_34835470的博客
03-17 592
通过为 cgs 提供函数句柄来求解线性方程组，该函数句柄计算 A*x 而不是系数矩阵 A。gallery 生成的 Wilkinson 测试矩阵之一是 21×21 三对角矩阵。预览该矩阵。A = gallery('wilk',21)A = 21×2110 1 0 0 0 0 0 0 0 0 0 0 0 ...
matlab 狄利克雷函数图像,Poisson image editing算法实现的Matlab代码解析
weixin_39609822的博客
03-17 2147
之前我发了数篇系列博文来仔细研究Poisson Image Editing算法，每次重新审视和深入，仿佛都能有更为深刻的认识和很大的收获。这应该算是我这个系列的完结篇，会用用Matlab代码一点一点的演示，原文作者到底是如何设计和实现他那个强大且影响深远的算法的。希望你在看本文之前务必参考一下文章来了解算法原理，本文将主要讲解编程实现的问题，对于前面讲过的内容，我不会深究。但我个人总体的感觉是，现...
MATLAB函数速查手册
落叶_19950624
03-25 1万+
内容简介 　　MATLAB是目前流行的理论与工程仿真软件之一。该软件自产生以来，就以其独有的特点和明显的优势吸引了各行各业的工作者。《MATLAB函数速查手册》较全面地介绍了MATLAB的函数，主要包括MATLAB操作基础、矩阵及其基本运算、与数值计算相关的基本函数、符号运算的函数、概率统计函数、绘图与图形处理函数、MATLAB程序设计相关函数、Simulink仿真工具函数、图形用户界面制作函数、...
matlab常用函数
good_learne的博客
09-07 1万+
一、软件操作函数 1）命令窗口函数： clc：清空命令窗口，使用向上箭头翻看命令。 open:打开文件，文本文件（*.doc），可执行文件（*.exe），图形文件（*.fig），超文本文件（*.html，*.htm），MATLAB数据库文件（*.mat），simulink模型文件*.mdl），MATLAB p文件（*.p），PDF文件（*,pdf），PPT文件（*.ppt），工程文件（*.prj...
matlab基本函数
freshman_kx的博客
04-05 7596
*matlab基本函数 1.生成矩阵 （1）直接法：a = [1,2,3;4,5,6;7,8,9] （2）冒号一维矩阵： a = 开始：步长：结束，步长为1可省略 b = 1:1:10; b = 1:10; （与上面等价） （3）函数生成： linspace(开始，结束，元素个数)，等差生成指定元素数的一维矩阵，省略个数则生成100个 c = linspace(0,10,5); （4）特殊矩阵： zeros函数:产生零矩阵 ones函数：产生全1矩阵 eye函数：产生对角线为1的矩阵 rand函数：
【整理】Matlab常用函数
qq_41633777的博客
05-13 1774
第一篇：Matlab软件函数 一、软件操作函数 1）命令窗口函数： clc：清空命令窗口，使用向上箭头翻看命令。 open:打开文件，文本文件（*.doc），可执行文件（*.exe），图形文件（*.fig），超文本文件（*.html，*.htm），MATLAB数据库文件（*.mat），simulink模型文...
Matlab从入门到精通（五）--常用函数
热门推荐
数据之道
11-23 2万+
1、傅立叶(Fourier)变换及拉普拉斯(Laplace)变换 在MATLAB中，进行傅立叶变换的函数是： fourier(f,x,t)：求函数f(x)的傅立叶像函数F(t)。 ifourier(F,t,x)：求傅立叶像函数F(t)的原函数f(x)。   在MATLAB中，进行拉普拉斯变换的函数是： laplace(fx,x,t)：求函数f(x)的拉普拉斯像函数F(t)。 ilap...
MATLAB 函数查询
MissXy_的博客
10-29 1万+
官方查询首页：http://cn.mathworks.com/help/matlab/functionlist.html#language-fundamentals MATLAB 函数 按字母顺序排列的列表按类别 语言基础知识 输入命令 ans 最近计算的答案 clc 清除命
matlab eigs 虚数,MATLAB函数速查手册
weixin_36019375的博客
03-16 1037
内容简介MATLAB是目前流行的理论与工程仿真软件之一。该软件自产生以来，就以其独有的特点和明显的优势吸引了各行各业的工作者。《MATLAB函数速查手册》较全面地介绍了MATLAB的函数，主要包括MATLAB操作基础、矩阵及其基本运算、与数值计算相关的基本函数、符号运算的函数、概率统计函数、绘图与图形处理函数、MATLAB程序设计相关函数、Simulink仿真工具函数、图形用户界面制作函数、信号处...
PhysConst:MATLAB 各种单位系统中基本物理常数的集合。-matlab开发
05-29
在MATLAB中使用PhysConst，首先需要下载并解压提供的`github_repo.zip`文件，然后按照标准方法将工具箱添加到MATLAB的工作路径中。一旦完成这些步骤，就可以通过简单的函数调用来获取所需常数。例如，若要获取普朗克...
matlab中实现IGBT仿真
03-15
### MATLAB中实现IGBT仿真的关键技术点 #### IGBT动态模型构建 在MATLAB中实现IGBT（绝缘栅双极晶体管）动态仿真，首先需要理解IGBT的基本结构及其动态特性。IGBT作为一种复合型功率器件，结合了MOSFET（金属...
Chapter10_Fig2_matlab_
10-03
MATLAB是一种广泛使用的编程环境，尤其在工程和科学计算中非常流行。它的语法简洁，适合进行数值计算、数据可视化以及编写复杂的算法。在这个场景下，MATLAB代码可能涉及到以下几个知识点： 1. **矩阵操作**：...
pingfanggen.rar_matlab数值求解_线性方程组
09-20
1. **`linsolve`函数**：这是MATLAB中解决一般线性方程组的标准方法，它能处理稠密或稀疏矩阵。当方程组的系数矩阵是方阵且可逆时，`linsolve`会使用高斯消元法或者其优化版本，如LU分解。 2. **`mldivide`（`\`)...
matlab开发-Jacobianfreenewtonkrylovjfnk方法
08-22
3. 迭代过程：在每次迭代中，使用Krylov子空间方法（如GMRES或CGS）求解线性化后的方程。 4. 更新解：根据求得的Krylov子空间解向量，更新当前解。 5. 检查收敛性：比较新的残差与旧的残差，如果满足收敛条件，则...
CODATA 2006：来自CODATA 2006的物理常数-matlab开发
05-30
**正文** ...使用这个类，科学家和工程师可以在MATLAB环境中进行精确的物理计算，而无需担心物理常数的准确性问题，从而提高研究和工程计算的可靠性。这体现了MATLAB的灵活性和对科学数据处理的强大支持。
matlab数学手册
01-11
- 使用`numel`函数计算矩阵中元素的数量。 - 示例：`X = numel(A);` **1.3 矩阵分解** - **1.3.1 Cholesky分解** - 使用`chol`函数进行Cholesky分解。 - 示例：`Y = chol(A);` - **1.3.2 LU分解** - 使用`lu...
Matlab数学手册
07-07
- 使用`cgs`函数。 - **1.4.8 共轭梯度的LSQR方法** - 使用`lsqr`函数。 - **1.4.9 广义最小残差法** - 使用`gmres`函数。 - **1.4.10 最小残差法解方程组** - 使用`minres`函数。 - **1.4.11 预处理共轭梯度...
Matlab中有单位吗
最新发布
03-30
Matlab中有单位，它支持SI、CGS、英制等多种单位。可以使用单位来进行数值计算和转换。Matlab中的单位可以通过符号表示，例如1米可以表示为1m，1秒可以表示为1s。Matlab还提供了一些预定义的单位，例如重力加速度可以表示为g0，光速可以表示为c。Matlab还提供了一些函数来进行单位转换，例如unitconvert函数可以将一个单位转换为另一个单位。
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值