求解线性系统

2024-07-14 07:37| 来源: 网络整理| 查看: 265

打开实时脚本

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

创建一个非对称三对角矩阵。预览该矩阵。

A = gallery('wilk',21) + diag(ones(20,1),1)A = 21×21 10 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 9 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 8 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 7 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 6 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 5 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 4 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 3 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 2 0 0 0 0 0 0 0 0 0 0 ⋮

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

表达式 A x 变为：

A x=[1020⋯⋯01920⋮01⋱20⋮010⋱⋱⋮0⋱1⋱0⋮⋱⋱⋱20⋯⋯0110][x1x2x3⋮⋮x21]=[10x1+2x2x1+9x2+2x3⋮⋮x19+9x20+2x21x20+10x21].

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

A x=[10x1+2x2x1+9x2+2x3⋮⋮x19+9x20+2x21x20+10x21]=[0x1x2⋮x20]+[10x19x2⋮9x2010x21]+2⋅[x2x3⋮x210]。

同样，AT x 的表达式变为：

AT x=[1010⋯⋯02910⋮02⋱10⋮020⋱⋱⋮0⋱1⋱0⋮⋱⋱⋱10⋯⋯0210][x1x2x3⋮⋮x21]=[10x1+x22x1+9x2+x3⋮⋮2x19+9x20+x212x20+10x21].

AT x=[10x1+x22x1+9x2+x3⋮⋮2x19+9x20+x212x20+10x21]=2⋅[0x1x2⋮x20]+[10x19x2⋮9x2010x21]+[x2x3⋮x210].

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

function y = afun(x,flag) if strcmp(flag,'notransp') % Compute A*x y = [0; x(1:20)] ... + [(10:-1:0)'; (1:10)'].*x ... + 2*[x(2:end); 0]; elseif strcmp(flag,'transp') % Compute A'*x y = 2*[0; x(1:20)] ... + [(10:-1:0)'; (1:10)'].*x ... + [x(2:end); 0]; end end

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

现在，通过为 qmr 提供用于计算 A*x 和 A'*x 的函数句柄，求解线性系统 Ax=b。使用容差 1e-6 和 25 次迭代。指定 b 为 A 的行总和，使得 x 的实际解是由 1 组成的向量。

b = full(sum(A,2)); tol = 1e-6; maxit = 25; x1 = qmr(@afun,b,tol,maxit)qmr converged at iteration 19 to a solution with relative residual 4.7e-07. x1 = 21×1 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 ⋮

局部函数

【本文地址】

求解线性系统

求解线性系统

今日新闻

推荐新闻