MATLAB中bicgstab函数用法

于 2024-08-20 08:35:11 发布
阅读量426 收藏 6
点赞数 5
分类专栏： Matlab 文章标签： matlab 算法 开发语言
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/jk_101/article/details/133163291
版权

目录

语法

说明

示例

线性系统的迭代解

使用指定了预条件子的 bicgstab

提供初始估计值

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

        bicgstab函数功能是求解线性系统 - 稳定双共轭梯度法。

语法

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

说明

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

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

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

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

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

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

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

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

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

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

示例

线性系统的迭代解

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

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

rng default
A = sprand(400,400,.5);
A = A'*A;
b = rand(400,1);

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

x = bicgstab(A,b);


bicgstab 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.12.

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

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

x = bicgstab(A,b,1e-4,100);


bicgstab stopped at iteration 100 without converging to the desired tolerance 0.0001
because the maximum number of iterations was reached.
The iterate returned (number 100) has relative residual 0.044.

        即使采用更宽松的容差和更多迭代，残差也并未改进多少。当迭代算法以这种方式停滞时，显然需要预条件子矩阵。

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

L = ichol(A);
x = bicgstab(A,b,1e-4,100,L');


bicgstab converged at iteration 30.5 to a solution with relative residual 5.3e-05.

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

使用指定了预条件子的 bicgstab

        检查使用指定了预条件子矩阵的 bicgstab 来求解线性系统的效果。

        加载 west0479，它是一个非对称的 479×479 实稀疏矩阵。

load west0479
A = west0479;

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

b = sum(A,2);

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

tol = 1e-12;
maxit = 20;

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

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

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

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

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

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

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

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

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

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

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

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

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')

如图所示：

提供初始估计值

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

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

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

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

maxit = 50;
x1 = bicgstab(A,b,[],maxit);
bicgstab converged at iteration 20.5 to a solution with relative residual 9.3e-07.
x0 = 0.99*e;
x2 = bicgstab(A,b,[],maxit,[],[],x0);
bicgstab converged at iteration 4 to a solution with relative residual 8.7e-07.

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

返回中间结果

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

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

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

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

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

        通过为 bicgstab 提供用来计算 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

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

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

b = ones(21,1);
tol = 1e-12;  
maxit = 50;
x1 = bicgstab(@afun,b,tol,maxit)
bicgstab converged at iteration 11.5 to a solution with relative residual 5.2e-13.
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

双共轭梯度稳定法

        双共轭梯度稳定 (BiCGSTAB) 算法是在双共轭梯度算法的基础上开发的，该算法通过使用重新启动的 GMRES 步骤来减轻 BiCG 步骤中的不规则收敛行为，从而在 BiCG 算法基础上进行改善。

        BiCGSTAB 中的 GMRES(1) 修复步骤仅利用一阶最小残差多项式，因此当这些低阶步骤不足时，该方法可能会停滞。更高阶多项式的使用导致了其他方法的开发，如 BiCGSTABL [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] van der Vorst, H.A., "BI-CGSTAB: A fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems," SIAM J. Sci. Stat. Comput., March 1992, Vol. 13, No. 2, pp. 631–644.

jk_101
关注
  • 5
    点赞
  • 踩
  • 6
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
Matlab在线性代数中的应用
05-19
SVD是线性代数中的另一个重要概念，Matlab中的`svd`函数可以进行奇异值分解。SVD在图像处理、数据分析等领域有广泛应用，如低秩恢复、矩阵分解等。 六、线性方程组的迭代求解 对于大型稀疏线性系统，直接求解方法...
C-N_matlab_CrankNicolson_
09-29
5. **求解线性系统**：使用MATLAB的内置函数如`sparse`创建稀疏矩阵，然后用`solve`或`bicgstab`等迭代方法求解线性系统 \(Au^{n+1}=b\)。 6. **时间迭代**：在每个时间步，重复4-5步，直到达到预设的终止时间。 7...
MATLAB中ilu函数用法
jk_101的博客
09-22 695
ilu函数的功能是对矩阵进行不完全 LU 分解。
MATLAB中bicg函数用法
jk_101的博客
08-19 339
对于 CG 算法所针对的对称正定方程组，BiCG 算法可提供相同的结果，但每次迭代的成本会增加一倍。如果 bicg 无法在达到最大迭代次数后收敛或出于任何原因暂停，则会显示一条包含相对残差 norm(b-A*x)/norm(b) 以及该方法停止时的迭代次数的诊断消息。绘制每个解的残差历史记录图，并添加一条表示指定容差的线。使用每行的总和作为 Ax=b 右侧的向量，使 x 的预期解是由 1 组成的向量。使用采用默认设置的 bicg 求解系数矩阵为方阵的线性系统，然后在求解过程中调整使用的容差和迭代次数。
matlab基本函数
freshman_kx的博客
04-05 7614
*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常用函数
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实用函数汇总
Newbeesmith的博客
05-02 3107
本文主要是一些实用函数，基础函数很少涉及，仅仅给出用法，也希望读者实际操作熟悉函数。 一、解方程组及矩阵变换 1.tril(A,n),triu(A,n) 求矩阵的上、下三角矩阵，A为原矩阵，n为主对角线平移量，默认为0，n>0向主对角线右上方平移 2.[L,U]=lu(A) 将矩阵A分解的上三角矩阵保存在U当中，将一个“心理学上的”下三角矩阵（例如一个下三角矩阵和置换矩阵的乘积）保存在L中，...
【整理】Matlab常用函数
qq_41633777的博客
05-13 1778
第一篇：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函数速查手册
落叶_19950624
03-25 1万+
内容简介 　　MATLAB是目前流行的理论与工程仿真软件之一。该软件自产生以来，就以其独有的特点和明显的优势吸引了各行各业的工作者。《MATLAB函数速查手册》较全面地介绍了MATLAB的函数，主要包括MATLAB操作基础、矩阵及其基本运算、与数值计算相关的基本函数、符号运算的函数、概率统计函数、绘图与图形处理函数、MATLAB程序设计相关函数、Simulink仿真工具函数、图形用户界面制作函数、...
MATLAB 函数查询
MissXy_的博客
10-29 1万+
官方查询首页：http://cn.mathworks.com/help/matlab/functionlist.html#language-fundamentals MATLAB 函数 按字母顺序排列的列表按类别 语言基础知识 输入命令 ans 最近计算的答案 clc 清除命
计算机数值分析课程matlab
11-21
2. 插值与拟合：学习使用多项式插值、样条插值和最小二乘拟合方法，了解插值误差分析，并在MATLAB中编写代码实现这些算法。 3. 微分方程数值解法：包括欧拉方法、龙格-库塔方法、隐式和显式差分法等，以及MATLAB的...
使用 MATLAB 编程学习 CFD 中的有限体积方法-matlab开发.zip
06-16
MATLAB 的迭代求解器，如 `bicgstab` 或 `gmres`，可用于求解隐式方法中的线性系统。对于显式方法，我们可以直接使用时间步长更新状态变量。 最后，结果后处理是 CFD 模拟不可或缺的一部分。MATLAB 提供了 `pcolor`...
SpTriSolve:稀疏三角求解-matlab开发
06-01
`SpTriSolve` 的主要优点在于，它专门优化了稀疏三角矩阵的求解过程，这在迭代方法中非常有用，例如用于预处理步骤，提高迭代求解的收敛速度。 迭代方法，如 GMRES、CG、BiCGStab 等，通常需要求解线性子问题，例如...
回归预测|基于灰狼GWO优化BP神经网络多输入多输出的数据回归预测Matlab程序GWO-BP 含预测新数据程序
2401_86241083的博客
08-19 806
GWO-BP（Grey Wolf Optimization and Backpropagation）模型结合了灰狼优化算法（GWO）和反向传播（BP）算法，旨在利用灰狼优化的全局搜索能力改进反向传播神经网络的训练过程。包括但不限于。
基于模糊神经网络的金融序列预测算法matlab仿真
熟悉python,matlab,C,C++,JAVA等
08-16 836
基于模糊神经网络的金融序列预测算法是一种结合了模糊逻辑和神经网络技术的先进预测方法，它适用于处理非线性、不确定性和模糊性的金融数据预测任务。传统的预测方法往往难以捕捉金融市场中的非线性关系和不确定性，而模糊神经网络因其独特的非线性映射能力和模糊逻辑处理能力，在处理这类问题上显示出优势。通过结合模糊逻辑的强大表达能力和神经网络的学习能力，这种算法能够捕捉到复杂的市场行为模式，从而为投资者提供更加准确的预测结果。一个典型的模糊规则可以表示为：“如果价格变化是高的并且MACD是正的，则预测的趋势是上涨”。
【数学建模】MATLAB快速入门
最新发布
2303_79015671的博客
08-19 648
本文将介绍MATLAB的一些基本用法如：矩阵的运算，做二维图，三维图等 更新中
MATLAB-基于Transformer模型的NASA锂电池剩余寿命预测
2401_86544683的博客
08-19 270
首先从NASA数据集中提取电池容量特征，然后基于B0005号电池数据训练，用B0006号电池数据测试预测
matlab bicgstab
08-17
BICGSTAB (BiConjugate Gradient STABilized) 是一种迭代解线性方程组的方法，常用于求解大型稀疏线性方程组。在Matlab中，bicgstab函数可以用来实现BICGSTAB算法。 你可以通过以下方式在Matlab中使用bicgstab函数来求解线性方程组： ```matlab % 定义系数矩阵A和右侧向量b A = ...; % 系数矩阵 b = ...; % 右侧向量 % 设置初始猜测解x0 x0 = ...; % 初始猜测解 % 调用bicgstab函数求解线性方程组 [x, flag, relres, iter, resvec] = bicgstab(A, b, [], [], [], [], x0); % 输出结果 disp('解:'); disp(x); disp('迭代次数:'); disp(iter); disp('相对残差:'); disp(relres); ``` 在上述代码中，A表示系数矩阵，b表示右侧向量，x0表示初始猜测解。bicgstab函数的输出包括解x、求解状态标志flag、相对残差relres、迭代次数iter和残差向量resvec。 请注意，具体使用时需要根据实际情况来定义系数矩阵A和右侧向量b。
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值