18.5专业的解算器

 
:x= bicg (A,b)
:x= bicg (A,b,tol)
:x= bicg (A,b,tol,maxit)
:x= bicg (A,b,tol,maxit,M)
:x= bicg (A,b,tol,maxit,M1,M2)
:x= bicg (A,b,tol,maxit,M, [],x0)
:x= bicg (A,b,tol,maxit,M1,M2,x0)
:x= bicg (A,b,tol,maxit,M, [],x0, …)
:x= bicg (A,b,tol,maxit,M1,M2,x0, …)
:[x,flag,relres,iter,resvec] = bicg (A,b, …)

求解线性方程组A * x = b 采用双共轭梯度迭代法。

输入参数为:

  • A是线性系统的矩阵,它必须是平方。A可以作为矩阵、函数句柄或内联函数传递Afcn使得Afcn (x, "notransp") = A * x 和Afcn (x, "transp") = A' * x 。的附加参数Afcn可能在之后通过x0.
  • b是右侧向量。它必须是行数与相同的列向量A.
  • tol是残余误差所需的相对误差范围,b - A * x 。如果规范(b - A * x)tol * norm (b)  如果tol则使用1e-6的误差范围。
  • maxit是允许的最大迭代次数;如果maxit则使用值20。
  • M1,M2是预处理器。预处理器M给定为M = M1 * M2二者都M1M2可以作为矩阵或函数句柄或内联函数传递g使得g (x, "notransp") = M1 \ x 或g (x, "notransp") = M2 \ x 和g (x, "transp") = M1' \ x 或g (x, "transp") = M2' \ x 如果M1省略或为空,则不应用预处理。预处理系统在理论上等效于应用bicg线性系统的方法inv (M1) * A * inv (M2) * y = inv (M1) * binv (M2') * A' * inv (M1') * z = inv (M2') * b然后设置x = inv (M2) * y.
  • x0是最初的猜测。如果x0被省略或为空,则函数集x0默认情况下为零向量。

以下任何参数x0被视为参数,并以适当的方式传递给任何函数(AfcnMfcn)或者已经给予bicg.

输出参数为:

  • x是的解的计算近似值A * x = b 。如果算法没有收敛,那么x是具有最小残差的迭代。
  • flag表示退出状态:
    • 0:算法收敛到规定的容差内。
    • 1:该算法没有收敛,并且达到了最大迭代次数。
    • 2:预处理矩阵是奇异的。
    • 3:算法停滞,即当前迭代之间的差的绝对值x而前面的小于eps * norm (x,2).
    • 4:算法无法继续,因为中间值太小或太大,无法进行可靠的计算。
  • relres是最终残差与其初始值的比率,以欧几里得范数测量。
  • iter是迭代x计算。
  • resvec是包含每次迭代时的残差的向量。执行的迭代总数从下式给出length (resvec) - 1.

考虑一个三对角矩阵的平凡问题

n = 20;
A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ...
    toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
              sparse (1, 2, 1, 1, n) * n / 2);
b = A * ones (n, 1);
restart = 5;
[M1, M2] = ilu (A);  # in this tridiag case, it corresponds to lu (A)
M = M1 * M2;
Afcn = @(x, string) strcmp (string, "notransp") * (A * x) + ...
                     strcmp (string, "transp") * (A' * x);
Mfcn = @(x, string) strcmp (string, "notransp") * (M \ x) + ...
                     strcmp (string, "transp") * (M' \ x);
M1fcn = @(x, string) strcmp (string, "notransp") * (M1 \ x) + ...
                     strcmp (string, "transp") * (M1' \ x);
M2fcn = @(x, string) strcmp (string, "notransp") * (M2 \ x) + ...
                     strcmp (string, "transp") * (M2' \ x);

示例1:的最简单用法bicg

x = bicg (A, b)

示例2: bicg具有一个计算A*xA'*x

x = bicg (Afcn, b, [], n)

示例3: bicg具有预处理矩阵M

x = bicg (A, b, 1e-6, n, M)

示例4: bicg具有作为预处理器的函数

x = bicg (Afcn, b, 1e-6, n, Mfcn)

示例5: bicg具有预处理矩阵M1M2

x = bicg (A, b, 1e-6, n, M1, M2)

示例6: bicg具有作为预处理器的函数

x = bicg (Afcn, b, 1e-6, n, M1fcn, M2fcn)

示例7: bicg将需要参数的函数作为输入

function y = Ap (A, x, string, z)
  ## compute A^z * x or (A^z)' * x
  y = x;
  if (strcmp (string, "notransp"))
    for i = 1:z
      y = A * y;
    endfor
  elseif (strcmp (string, "transp"))
    for i = 1:z
      y = A' * y;
    endfor
  endif
endfunction

Apfcn = @(x, string, p) Ap (A, x, string, p);
x = bicg (Apfcn, b, [], [], [], [], [], 2);

参考

Y.Saad,Iterative Methods for Sparse Linear Systems,2003年第二版,SIAM。

详见: bicgstab,cgs,gmres,pcg,qmr,tfqmr.

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

解决A x = b使用稳定的双共轭梯度迭代方法。

输入参数为:

以下参数x0被视为参数,并以适当的方式传递给任何函数(AM)传递给bicstab.

输出参数为:

让我们考虑一个三对角矩阵的平凡问题

n = 20;
A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n))  + ...
    toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
    sparse (1, 2, 1, 1, n) * n / 2);
b = A * ones (n, 1);
restart = 5;
[M1, M2] = ilu (A); # in this tridiag case, it corresponds to lu (A)
M = M1 * M2;
Afcn = @(x) A * x;
Mfcn = @(x) M \ x;
M1fcn = @(x) M1 \ x;
M2fcn = @(x) M2 \ x;

示例1:的最简单用法bicgstab

x = bicgstab (A, b, [], n)

示例2: bicgstab具有一个计算A * x

x = bicgstab (Afcn, b, [], n)

示例3: bicgstab具有预处理矩阵M

x = bicgstab (A, b, [], 1e-06, n, M)

示例4: bicgstab具有作为预处理器的函数

x = bicgstab (Afcn, b, 1e-6, n, Mfcn)

示例5: bicgstab具有预处理矩阵M1M2

x = bicgstab (A, b, [], 1e-6, n, M1, M2)

示例6: bicgstab具有作为预处理器的函数

x = bicgstab (Afcn, b, 1e-6, n, M1fcn, M2fcn)

示例7: bicgstab将一个需要参数的函数作为输入

function y = Ap (A, x, z) # compute A^z * x
   y = x;
   for i = 1:z
     y = A * y;
   endfor
 endfunction
Apfcn = @(x, string, p) Ap (A, x, string, p);
x = bicgstab (Apfcn, b, [], [], [], [], [], 2);

示例8:明确的例子表明bicgstab使用正确的预处理器

[M1, M2] = ilu (A + 0.1 * eye (n)); # factorization of A perturbed
M = M1 * M2;

## reference solution computed by bicgstab after one iteration
[x_ref, fl] = bicgstab (A, b, [], 1, M)

## right preconditioning
[y, fl] = bicgstab (A / M, b, [], 1)
x = M \ y # compare x and x_ref

参考

Y.Saad,Iterative Methods for Sparse Linear Systems,2003年第二版,SIAM

详见: bicg,cgs,gmres,pcg,qmr,tfqmr.

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

解决A x = b这里的A是一个正方阵,使用共轭梯度平方法。

输入参数为:

以下参数x0被视为参数,并以适当的方式传递给任何函数(AP)传递给cgs.

输出参数为:

让我们考虑一个三对角矩阵的平凡问题

n = 20;
A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n))  + ...
    toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
    sparse (1, 2, 1, 1, n) * n / 2);
b = A * ones (n, 1);
restart = 5;
[M1, M2] = ilu (A); # in this tridiag case it corresponds to chol (A)'
M = M1 * M2;
Afcn = @(x) A * x;
Mfcn = @(x) M \ x;
M1fcn = @(x) M1 \ x;
M2fcn = @(x) M2 \ x;

示例1:的最简单用法cgs

x = cgs (A, b, [], n)

示例2: cgs具有一个计算A * x

x = cgs (Afcn, b, [], n)

示例3: cgs具有预处理矩阵M

x = cgs (A, b, [], 1e-06, n, M)

示例4: cgs具有作为预处理器的函数

x = cgs (Afcn, b, 1e-6, n, Mfcn)

示例5: cgs具有预处理矩阵M1M2

x = cgs (A, b, [], 1e-6, n, M1, M2)

示例6: cgs具有作为预处理器的函数

x = cgs (Afcn, b, 1e-6, n, M1fcn, M2fcn)

示例7: cgs将需要参数的函数作为输入

function y = Ap (A, x, z) # compute A^z * x
   y = x;
   for i = 1:z
     y = A * y;
   endfor
 endfunction
Apfcn = @(x, string, p) Ap (A, x, string, p);
x = cgs (Apfcn, b, [], [], [], [], [], 2);

示例8:明确的例子表明cgs使用正确的预处理器

[M1, M2] = ilu (A + 0.3 * eye (n)); # factorization of A perturbed
M = M1 * M2;

## reference solution computed by cgs after one iteration
[x_ref, fl] = cgs (A, b, [], 1, M)

## right preconditioning
[y, fl] = cgs (A / M, b, [], 1)
x = M \ y # compare x and x_ref

参考文献:

Y.Saad,Iterative Methods for Sparse Linear Systems,2003年第二版,SIAM

详见: pcg,bicgstab,bicg,gmres,qmr,tfqmr.

广告
 
:x= gmres (A,b,restart,tol,maxit,M1,M2,x0, …)
:x= gmres (A,b,restart,tol,maxit,M, [],x0, …)
:[x,flag,relres,iter,resvec] = gmres (A,b, …)

解决A x = b使用带重启的预处理GMRES迭代方法,也称为PGMRES(重启)。

输入参数为:

以下参数x0被视为参数,并以适当的方式传递给任何函数(AMM1M2)传递到gmres.

输出为:

让我们考虑一个三对角矩阵的平凡问题

n = 20;
A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n))  + ...
    toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
    sparse (1, 2, 1, 1, n) * n / 2);
b = A * ones (n, 1);
restart = 5;
[M1, M2] = ilu (A); # in this tridiag case, it corresponds to lu (A)
M = M1 * M2;
Afcn = @(x) A * x;
Mfcn = @(x) M \ x;
M1fcn = @(x) M1 \ x;
M2fcn = @(x) M2 \ x;

示例1:的最简单用法gmres

x = gmres (A, b, [], [], n)

示例2: gmres具有一个计算A * x

x = gmres (Afcn, b, [], [], n)

示例3:的用法gmres随着重新启动

x = gmres (A, b, restart);

示例4: gmres具有预处理矩阵M有无重新启动

x = gmres (A, b, [], 1e-06, n, M)
x = gmres (A, b, restart, 1e-06, n, M)

示例5: gmres具有作为预处理器的函数

x = gmres (Afcn, b, [], 1e-6, n, Mfcn)

示例6: gmres具有预处理矩阵M1M2

x = gmres (A, b, [], 1e-6, n, M1, M2)

示例7: gmres具有作为预处理器的函数

x = gmres (Afcn, b, 1e-6, n, M1fcn, M2fcn)

示例8: gmres将需要参数的函数作为输入

  function y = Ap (A, x, p) # compute A^p * x
     y = x;
     for i = 1:p
       y = A * y;
     endfor
  endfunction
Apfcn = @(x, p) Ap (A, x, p);
x = gmres (Apfcn, b, [], [], [], [], [], [], 2);

示例9:明确的例子表明gmres使用aleft预处理器

[M1, M2] = ilu (A + 0.1 * eye (n)); # factorization of A perturbed
M = M1 * M2;

## reference solution computed by gmres after two iterations
[x_ref, fl] = gmres (A, b, [], [], 1, M)

## left preconditioning
[x, fl] = gmres (M \ A, M \ b, [], [], 1)
x # compare x and x_ref

参考

Y.Saad,Iterative Methods for Sparse Linear Systems,2003年第二版,SIAM

详见: bicg,bicgstab,cgs,pcg,pcr,qmr,tfqmr.

广告
 
:x= qmr (A,b,rtol,maxit,M1,M2,x0)
:x= qmr (A,b,rtol,maxit,P)
:[x,flag,relres,iter,resvec] = qmr (A,b, …)

解决A x = b使用准最小残差迭代方法(无需前瞻)。

A可以作为矩阵或函数句柄或内联函数传递f使得f(x, "notransp") = A*xf(x, "transp") = A'*x.

预处理器P给定为P = M1 * M2二者都M1M2可以作为矩阵或函数句柄或内联函数传递g使得g(x, "notransp") = M1 \ xg(x, "notransp") = M2 \ xg(x, "transp") = M1' \ xg(x, "transp") = M2' \ x.

如果使用多个输出参数调用

参考文献:

  1. R.Freund和N.Nachtigal,QMR: a quasi-minimal residual method for non-Hermitian linear systems,Numeriche Mathematik,1991,60,第315–339页。
  2. R.Barrett、M.Berry、T.Chan、J.Demmel、J.Donato、J.栋加拉、V.Eijkhour、R.Pozo、C.Romine和H.van der Vorst,Templates for the solution of linear systems: Building blocks for iterative methods,SIAM,第2版,1994年。

详见: bicg,bicgstab,cgs,gmres,pcg.

广告
 
:x= tfqmr (A,b,tol,maxit,M1,M2,x0, …)
:x= tfqmr (A,b,tol,maxit,M, [],x0, …)
:[x,flag,relres,iter,resvec] = tfqmr (A,b, …)

解决A x = b使用基于cgs的Transpose Tree qmr方法。

输入参数为:

以下参数x0被视为参数,并以适当的方式传递给任何函数(AM)传递给tfqmr.

输出参数为:

让我们考虑一个三对角矩阵的平凡问题

n = 20;
A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n))  + ...
    toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ...
    sparse (1, 2, 1, 1, n) * n / 2);
b = A * ones (n, 1);
restart = 5;
[M1, M2] = ilu (A); # in this tridiag case it corresponds to chol (A)'
M = M1 * M2;
Afcn = @(x) A * x;
Mfcn = @(x) M \ x;
M1fcn = @(x) M1 \ x;
M2fcn = @(x) M2 \ x;

示例1:的最简单用法tfqmr

x = tfqmr (A, b, [], n)

示例2: tfqmr具有一个计算A * x

x = tfqmr (Afcn, b, [], n)

示例3: tfqmr具有预处理矩阵M

x = tfqmr (A, b, [], 1e-06, n, M)

示例4: tfqmr具有作为预处理器的函数

x = tfqmr (Afcn, b, 1e-6, n, Mfcn)

示例5: tfqmr具有预处理矩阵M1M2

x = tfqmr (A, b, [], 1e-6, n, M1, M2)

示例6: tfmqr具有作为预处理器的函数

x = tfqmr (Afcn, b, 1e-6, n, M1fcn, M2fcn)

示例7: tfqmr将需要参数的函数作为输入

function y = Ap (A, x, z) # compute A^z * x
   y = x;
   for i = 1:z
     y = A * y;
   endfor
 endfunction
Apfcn = @(x, string, p) Ap (A, x, string, p);
x = tfqmr (Apfcn, b, [], [], [], [], [], 2);

示例8:明确的例子表明tfqmr使用正确的预处理器

[M1, M2] = ilu (A + 0.3 * eye (n)); # factorization of A perturbed
M = M1 * M2;

## reference solution computed by tfqmr after one iteration
[x_ref, fl] = tfqmr (A, b, [], 1, M)

## right preconditioning
[y, fl] = tfqmr (A / M, b, [], 1)
x = M \ y # compare x and x_ref

参考

Y.Saad,Iterative Methods for Sparse Linear Systems,2003年第二版,SIAM

详见: bicg,bicgstab,cgs,gmres,pcg,qmr,pcr.

广告

版权所有 © 2024-2025 Octave中文网

ICP备案/许可证号:黑ICP备2024030411号-2