精通MATLAB科学计算（第3版）(王正林)09-3r.pdf

第 章 线 性 方 程 组 求 解 在 自 然 科 学 和 工 程 技 术 中，很 多 问 题 可 以 归 结 为 求 解 线 性 方 程 组。采 用 MATLAB,不 仅 可 以 利 用 其 提 供 的 相 关 函 数 直 接 解 决 一 些 简 单 的 线 性 方 程 组，而 且 可 以 通 过 简 洁 的 编 程 来 解 决 一 些 复 杂 的 线 性 方 程 组。通 过 本 章 的 介 绍，读 者 既 能 应 用 MATLAB中 相 应 的 函 数 求 解 线 性 方 程 组，又 能 通 过 编 程，灵 活 使 用 迭 代 法 和 其 他 的 特 殊 解 法 来 求 解 线 性 方 程 组。9.1 求 逆 法 MATLAB中 求 解 线 性 方 程 最 直 接 的 方 法 是 矩 阵 求 逆 法，它 适 用 于 系 数 矩 阵 的 数 据 无 规 律 且 系 数 矩 阵 的 阶 数 比 较 小 的 情 况。如 果 系 数 矩 阵 的 阶 数 太 大 的 话，系 数 矩 阵 求 逆 就 需 要 花 费 很 长 时 间。在 MATLAB中 用 求 逆 法 求 解，线 性 方 程 可 以 直 接 用 左 除 法“”求 解，也 可 以 用 求 逆 函 数 inv()求 解。【例 9-1】左 除 法 和 求 逆 法 求 解 线 性 方 程 组 实 例。用 左 除 法 和 求 逆 法 求 解 如 下 线 性 方 程 组，X+2%2+3x3=1-x+3工 2+7工 3=49x i+3x3=7解：用 左 除 法 求 解，在 MATLAB命 令 窗 口 中 输 入 求 解 程 序：A=l 2 3;-1 3 7;9 0 3 b=l 4 7 z;x=A b输 出 计 算 结 果 为：X=0.3 3 3 3-1.6 6 6 71.3 3 3 3用 逆 矩 阵 法 求 解，在 MATLAB命 令 窗 口 中 输 入 求 解 程 序：第 9 章 线 性 方 程 组 求 解 A=l 2 3;-1 3 7;9 0 3;b=l 4 7 z;x=in v(A)*b输 出 计 算 结 果 为：X=0.3 3 3 3-1.6 6 6 71.3 3 3 3由 计 算 结 果 可 知,方 程 组 的 解 为 的 4 2户 3=03333,-1.6667,1.3333。9.2 分 解 法 矩 阵 分 解 法 是 指 根 据 一 定 的 原 理 用 某 种 算 法 将 系 数 矩 阵 分 解 成 若 干 个 矩 阵 的 代 数 运 算，常 用 的 分 解 是 乘 积 分 解，而 分 解 后 的 新 矩 阵 一 般 是 某 种 特 殊 矩 阵。常 用 的 分 解 有 L U分 解、Q R分 解、Cholesky分 解、Schur分 解、Hessenberg分 解 和 奇 异 分 解 等。9.2.1 LU分 解 法 矩 阵 的 L U分 解 也 叫 三 角 分 解，还 有 的 书 称 之 为 Doolittle分 解，它 是 将 矩 阵 分 解 为 一 个 单 位 下 三 角 矩 阵 与 上 三 角 矩 阵 的 乘 积。只 要 矩 阵 非 奇 异，这 种 分 解 总 是 可 以 进 行 的。MATLAB提 供 了 lu函 数 来 求 矩 阵 A 的 L U分 解，常 用 格 式 为：也,S=lu(Z)：产 生 一 个 上 三 角 矩 阵 U 和 一 个 变 换 形 式 的 下 三 角 矩 阵 Z(进 行 了 行 交 换)，使 得 A=LU；,(/,P=lu(/1)：产 生 一 个 上 三 角 矩 阵 U 和 一 个 下 三 角 矩 阵 乙 及 置 换 矩 阵 P,使 得 PA=LU；分 解 以 后，方 程 组=的 解 可 写 为 x=LA(ZA)或 x=(A(ZAP*)。【例 9-21 L U分 解 法 求 解 线 性 方 程 组 实 例。用 L U分 解 法 求 解 以 下 线 性 方 程 组，1.5xi+3x2-0.83+414=42xi+9%3+10%4=0 lx+4.8x2 0.6x3+=11 4%|+12.3%2 4M+5X 4=2 解：用 L U分 解 法 求 解，在 MATLAB命 令 窗 口 中 输 入 求 解 程 序：A=1.5 3-0.8 4;2 0 9 1 0;-7 4.8-0.6 1;1 4 1 2.3-4 5;b=4 0 1-2;Lz U=lu(A);1 187精 通 M ATLAB科 学 计 算(第 2 忡-x=U(L b)输 出 计 算 结 果 为：x=-0.3 7 2 1-0.8 2 9 1-1.5 3 3 61.4 5 4 7由 计 算 结 果 可 知,方 程 组 的 解 为 即 X2,X3,X4=-0.3721,-0.8291,-1.5336,1.4547。9.2.2 QR分 解 法 矩 阵 的 Q R 分 解 就 是 把 矩 阵 分 解 为 一 个 正 交 矩 阵 和 一 个 上 三 角 矩 阵 的 乘 积。MATLAB中 对 矩 阵 A 进 行 Q R 分 解 的 函 数 调 用 格 式 如 下：Q,K=qr(/)：产 生 一 个 正 交 矩 阵。和 一 个 上 三 角 矩 阵 A,使 得；Q,K,E=qr(/)：产 生 一 个 正 交 矩 阵 Q、一 个 上 三 角 矩 阵 R 以 及 一 个 置 换 矩 阵 E,使 得 AE=QR；实 现 分 解 以 后，方 程 组 4r=6 的 解 可 写 为 x=K(Q功)或 x=E(K(QS)。【例 9-31 Q R 分 解 法 求 解 线 性 方 程 组 实 例。用 Q R 分 解 法 求 解 以 下 线 性 方 程 组，x+0.5x2+0.3333%3+0.25%4=10.5%1+0.3333%2+0.25x3+0.2%4=20.3333x1+0.25x2+0.2%3+0.1667x4=20.25xi+0.2必+0.1667x3+0.1429x4=1。解：用 Q R 分 解 法 求 解，在 M A T L A B 命 令 窗 口 中 输 入 求 解 程 序：A=1.0 0 0 00.5 0 0 00.3 3 3 30.2 5 0 00.5 0 0 00.3 3 3 30.2 5 0 00.2 0 0 00.3 3 3 30.2 5 0 00.2 0 0 00.1 6 6 70.2 5 0 0;0.2 0 0 0;0.1 6 6 7;0.1 4 2 9;b=l 2 2 1,;Q,R=q r(A)x=R(Q b)输 出。、A 的 计 算 结 果 为：Q=-0.8 3 8 1 0.5 2 2 6-0.1 5 4 0-0.0 2 6 3-0.4 1 9 1-0.4 4 1 7 0.7 2 7 8 0.3 1 5 7-0.2 7 9 4-0.5 2 8 8-0.1 3 9 5-0.7 8 9 2-0.2 0 9 5-0.5 0 2 1-0.6 5 3 6 0.5 2 6 1-1.1 9 3 2-0.6 7 0 5-0.4 7 4 9-0.3 6 9 80-0.1 1 8 5-0.1 2 5 7-0.1 1 7 5188 第 9 章 线 性 方 程 组 求 解 0 0-0.0062-0.00960 0 0 0.0002输 出 解 的 计 算 结 果 为：x=1.0e+003*0.1160-1.44003.6000-2.3800由 计 算 结 果 可 知,方 程 组 的 解 为 的 广 2,孙 孙=口 的,-1440,3600,-2380。9.2.3 Cholesky 分 解 法 当 系 数 矩 阵 N 正 定 且 对 称 时，它 能 分 解 为 以 下 的 形 式:其 中 K 为 上 三 角 矩 阵，肝 为 R 的 转 置，是 下 三 角 矩 阵，这 种 分 解 称 为 Cholesky分 解。MATLAB中 与 Cholesky分 解 有 关 的 函 数 如 下：R=chol(Z):对 进 行 Cholesky 分 解，使 得 A=RXR;K,p=chol(4):对 A 进 行 Cholesky 分 解，使 得 A=RrRa如 果 4 对 称 正 定，则 p=0；否 则，p 为 一 正 整 数。如 果 4 满 秩，/?是 为 2-1 阶 的 上 三 角 矩 阵，且 有 相/?=4 1:仍-1),l:(p-l)；实 现 分 解 后，方 程 组 的 解 写 成 上 大(如 财 的 形 式。【例 9-4 Cholesky分 解 法 求 解 线 性 方 程 组 实 例。用 Cholesky分 解 法 求 解 以 下 线 性 方 程 组，9xi-36切+3013=1*-36xi+1922 一 180工 3=130%1-180%2+180%3=1 0解：用 Cholesky分 解 法 求 解，在 MATLAB命 令 窗 口 中 输 入 求 解 程 序：A=9-36 30;-36 192-180;30-180 180);b=ones(3,1);R=chol(A)x=R(Rzb)输 出/?的 计 算 结 果 为：R=3.0000-12.0000 10.00000 6.9282-8.66030 0 2.2361输 出 解 的 计 算 结 果 为：189精 通 MATLAB科 学 计 算(第 2 版 i-X=1.8 3 3 31.0 8 3 30.7 8 3 3由 计 算 结 果 可 知，方 程 组 的 解 为 莺/2,必=1.8333,1.0833,0.7833。起 合 在 M A TLA B中，当 N 正 定 但 不 对 称 时，用 K=chol(N)也 能 得 出 一 个 矩 阵，MATLAB并 不 报 错，但 是 只 要 验 证 一 下 Rt*R就 会 看 出 R*R 并 不 等 于 A,因 此 用 公 式 x=R(*M)就 会 得 出 错 误 的 结 果。190 第 9 章 线 性 方 程 组 求 解 9.2.4 其 他 分 解 法 MATLAB中 还 提 供 了 矩 阵 的 奇 异 值 分 解、Hessenberg分 解 和 Schur分 解 等 函 数，可 用 来 求 解 线 性 方 程 组，如 表 9-1所 示。表 9-1 其 他 的 矩 阵 分 解 函 数 函 数 及 常 用 语 法 说 明 U SM=svd(j)奇 异 值 分 解，使 得 P.H=hess()Hessenberg分 解.使 得 A=P*H*/，其 中 P 为 酉 矩 阵 U,T=schur(/l)Schui分 解，使 得 力。，其 中 U 为 正 交 矩 阵，T为 Schur矩 阵【例 9-5】奇 异 值 分 解 法 求 解 线 性 方 程 组 实 例。用 奇 异 值 分 解 法 求 解 以 下 线 性 方 程 组，X)+0.5x2+0.3333%3+0.25x4+0.2x5=10.5%1+0.3333x2+0.25%3+0.2x4+0.1667x5=0 UzSrV=svd(A)b=l 0 1 0 1 zr x=V*in v(S)*UI*b输 出 分 解 的 I/、S、H 的 计 算 结 果 为：U=-0.7 6 7 9 0.6 0 1 9-0.2 1 4 2 0.0 4 7 2 0.0 0 6 2-0.4 4 5 8-0.2 7 5 9 0.7 2 4 1-0.4 3 2 7-0.1 1 6 7-0.3 2 1 6-0.4 2 4 9 0.1 2 0 5 0.6 6 7 4 0.5 0 6 2-0.2 5 3 4-0.4 4 3 9-0.3 0 9 6 0.2 3 3 0-0.7 6 7 2-0.2 0 9 8-0.4 2 9 0-0.5 6 5 2-0.5 5 7 6 0.3 7 6 2S=1.5 6 7 1 0 0 0 00 0.2 0 8 5 0 0 00 0 0.0 1 1 4 0 00 0 0 0.00 0 3 00 0 0 0 0.,0000V=-0.7 6 7 9 0.6 0 1 9-0.2 1 4 2 0.0 4 7 2 0.0 0 6 2 191精 通 MATLAB科 学 计 算（第 2 蚪-0.4458-0.3216-0.2534-0.2098-0.2759-0.4249-0.4439-0.42900.72410.1205-0.3096-0.5652-0.43270.66740.2330-0.5576-0.11670.5062-0.76720.3762输 出 的 计 算 结 果 为：x=1.0e+004*0.0865-1.24674.4659-5.84042.5426由 计 算 结 果 可 知，方 程 组 的 解 为，x1,X2,x3,X4,x5=104*0.0865,-1.2467,4.4659,-5.8404,2.5426。选 W 奇 异 值 分 解 很 有 用，将 系 数 矩 阵 进 行 奇 异 值 分 解 以 后，A 的 奇 异 值 都 在 S的 对 角 线 上，这 样 就 可 以 粗 略 估 计 N 的 条 件 数，看 看 N 是 否 是 病 态 矩 阵，以 此 决 定 相 应 的 求 解 方 法。如 果 系 数 矩 阵 对 称，那 么 奇 异 值 分 解 后 的。和/都 是 正 交 矩 阵，而 S 是 对 角 矩 阵，这 样 求 解 方 程 就 十 分 方 便。【例 9-6 Hessenberg分 解 法 求 解 线 性 方 程 组 实 例。用 Hessenberg分 解 法 求 解 以 下 线 性 方 程 组，闲+工 2+工 3+工 4=1X+2X2+3汹+414=4X+3X 2+6x3+1014=7X+4X 2+1013+204=-2 o解：用 Hessenberg分 解 法 求 解，在 MATLAB命 令 窗 口 中 输 入 求 解 程 序：A=1 1 1 1;1 2 3 4;1 3 6 10;1 4 10 20;P,H=hess(A)b=l 4 7-21;x=P*inv(H)*P1*b输 出 分 解 的 P、的 计 算 结 果 为：P=0.7754-0.6247-0.0925 0-0.6092-0.7015-0.3698 0192 第 9 章 线 性 方 程 组 求 解 0.166200.34310-0.92450 1.0,00000.2147-0.2654 0 0-0.2654 1.0844 1.0802 00 1.0802 7.7009-10.81670 0-10.8167 20.0000输 出 的 计 算 结 果 为：10.0000-33.000036.0000-12.0000由 计 算 结 果 可 知,方 程 组 的 解 为,x i,x2,X 3=10,-3 3,3 6,-1 2 o起 合 Hessenberg矩 阵 指 的 是 第 一 子 对 角 线 以 下 的 元 素 为 0 的 矩 阵。Hessenberg分 解 法 分 解 后 产 生 正 交 矩 阵 和 Hessenberg矩 阵，正 交 矩 阵 的 逆 就 是 其 转 置；如 果 系 数 矩 阵 对 称 的 话，产 生 的 Hessenberg矩 阵 就 是 三 对 角 矩 阵，这 种 矩 阵 的 求 逆 速 度 也 很 快，因 此 用 这 种 方 法 求 解 线 性 方 程 组 也 很 快。【例 9-7】S chur分 解 法 求 解 线 性 方 程 组 实 例。S chur分 解 法 求 解 以 下 线 性 方 程 组，汨+工 2+工 3+工 4=1X+2X2+3冷+414=4X 1+3x2+6x3+10 x4=7X+42+10工 3+20%4=-2 o解：用 S chur分 解 法 求 解，在 M A T L A B命 令 窗 口 中 输 入 求 解 程 序:A=1 1 1 1；1 2 3 4;1 3 6 10;1 4 10 20;b=l 4 7-2 UzT=schur(A)x=U*inv(T)*Uf*b输 出 分 解 的 U、7 的 计 算 结 果 为：U=0.3087-0.7873 0.5304 0.0602-0.7231 0.1632 0.6403 0.20120.5946 0.5321 0.3918 0.4581-0.1684-0.2654-0.3939 0.8638T=0.0380 0 0 00 0.4538 0 01 1939.3精 通 M ATLAB科 学 计 算（第 2 版 i-0 0 2.2034 00 0 0 26.3047输 出 的 计 算 结 果 为：x=10.0000-33.000036.0000-12.0000由 计 算 结 果 可 知,方 程 组 的 解 为 XI,X2,X3,X4=10,-33,36,-12。通 多 Schur矩 阵 是 上 三 角 矩 阵，且 其 对 角 元 素 为 被 分 解 矩 阵 的 特 征 值。Schur分 解 法 在 分 解 后 产 生 正 交 矩 阵 和 Schur矩 阵，如 果 系 数 矩 阵 对 称 的 话，产 生 的 Schur矩 阵 就 是 对 角 阵，显 然 对 角 阵 的 逆 矩 阵 非 常 容 易 得 到，因 此 用 这 种 方 法 求 解 线 性 方 程 组 也 很 快。迭 代 法 在 实 际 应 用 中（例 如 在 运 筹 学、图 论 等 领 域 中）,往 往 会 出 现 这 样 一 种 情 况：系 数 矩 阵 阶 数 很 高（例 如 IC T）,但 系 数 矩 阵 含 零 元 素 相 对 较 多，如 果 用 前 面 的 几 种 分 解 法 去 求 解 的 话，非 零 元 素 反 而 增 多 了，这 就 有 点 得 不 偿 失 了。本 节 介 绍 另 一 类 常 用 的 求 解 方 法 迭 代 法。迭 代 法 是 将 求 一 组 解 转 换 为 求 一 个 近 似 解 序 列 的 过 程，并 用 最 终 的 近 似 解 来 逼 近 真 实 解。迭 代 法 需 要 考 虑 以 下 3 个 重 要 的 问 题。（1）迭 代 的 初 始 值 考 虑 初 始 值 的 选 取 是 否 有 范 围 限 制，不 同 的 初 始 值 对 最 终 的 迭 代 结 果 是 否 有 影 响。（2）迭 代 算 法 考 虑 怎 样 由 当 前 的 迭 代 结 果 得 出 下 次 迭 代 的 初 始 值；由 于 不 同 的 算 法 会 带 来 不 同 的 误 差，所 以 应 考 虑 算 法 使 误 差 放 大 还 是 减 少。（3）迭 代 的 收 敛 性 考 虑 迭 代 是 否 收 敛，收 敛 的 过 程 是 快 还 是 慢。以 上 问 题 的 具 体 分 析 请 参 考 数 值 分 析 的 教 材，下 面 讲 述 采 用 MATLAB编 程 实 现 的 几 种 常 见 迭 代 法。9.3.1 逐 次 逼 近 法 对 于 n 阶 线 性 方 程 组，Ax=b的 系 数 矩 阵 N（假 设 N 为 非 奇 异）进 行 如 下 分 解 A=Q-C194 第 9 章 线 性 方 程 组 求 解 其 中。为 非 奇 异。则 方 程 组 可 变 换 为 x=Bx+r其 中 B=Q C,r=Q：b,那 么 迭 代 过 程 可 以 写 成 xk+l=Bxk+r这 种 迭 代 过 程 称 为 逐 次 逼 近 法。取 不 同 的。和 C,就 得 到 不 同 的 迭 代 算 法。通 常 情 况 下，逐 次 逼 近 法 收 敛 的 充 分 条 件 是 迭 代 矩 阵 的 谱 半 径 小 于 lo9.3.2 里 查 森 迭 代 法 里 查 森 迭 代 法 可 以 说 是 最 简 单 的 迭 代 法 了，它 采 用 的 迭 代 公 式 如 下：xk+i=(I-A)xk+b在 MATLAB中 编 程 实 现 的 里 查 森 迭 代 法 函 数 为：richasono功 能：用 里 查 森 迭 代 法 求 线 性 方 程 组 的 解。调 用 格 式：x,n=richason(,b,xO,eps,M)其 中，A 为 线 性 方 程 组 的 系 数 矩 阵；分 为 线 性 方 程 组 中 的 常 数 向 量；xO为 迭 代 初 始 向 量；e p s为 解 的 精 度 控 制(此 参 数 可 选)；M 为 迭 代 步 数 控 制(此 参 数 可 选)；x 为 线 性 方 程 组 的 解；为 求 出 所 需 精 度 的 解 实 际 的 迭 代 步 数。理 查 森 迭 代 法 的 MATLAB程 序 代 码 如 下：function x,n=richason(A,b,xO,eps,M)%采 用 里 查 森 迭 代 法 求 线 性 方 程 组 Ax=b的 解 W线 性 方 程 组 的 系 数 矩 阵：A告 线 性 方 程 组 中 的 常 数 向 量：b%迭 代 初 始 向 量：xO3解 的 精 度 控 制：eps为 迭 代 步 数 控 制：M年 线 性 方 程 组 的 解：x若 求 出 所 需 精 度 的 解 实 际 的 迭 代 步 数：nif(nargin=3)eps=1.Oe-6;%eps表 示 迭 代 精 度 1 195精 通 MATLAB科 学 计 算(第 2 忡-M=200;*M 表 示 迭 代 步 数 的 限 制 值 elseif(nargin=4)M=200;endI=eye(size(A);xl=x0;x=(I-A)*x0+b;n=l;告 迭 代 过 程 while(norm(x-xl)eps)xl=x;x=(I-A)*xl+b;n=n+1;会 n 为 最 终 求 出 解 时 的 迭 代 步 数 if(n=M)disp(Warning:迭 代 次 数 太 多，可 能 不 收 敛！,);return;endend【例 9 8】理 查 森 迭 代 法 求 解 线 性 方 程 组 实 例。用 理 查 森 迭 代 法 求 解 以 下 线 性 方 程 组，其 中 取 初 始 值 为 0,0,0。1.0170汨 一 0.0092M-0.0095x3=1,-0.0092X+0.9903x2+0.0136x3=0-0.0095xi+0.0136x2+0.9898x3=1 0解：用 理 查 森 迭 代 法 求 解，在 MATLAB命 令 窗 口 中 输 入 求 解 程 序：A=1.0170-0.0092 0.0095;-0.0092 0.9903 0.0136;0.0095 0.0136 0.9898;b=l 0 1;x0=0 0 0/_;x,n=richason(A,b,xO)输 出 的 计 算 结 果 为：X=0.9738-0.00471.0010输 出 的 迭 代 次 数 为：n=g身 可 见，经 过 5 步 迭 代，理 查 森 迭 代 法 求 出 了 方 程 组 的 解 为：196 第 9 章 线 性 方 程 组 求 解 阳,xz,X3=0.9738；6.0047,1.0010 l巧，巧，演=0.9 7 3 8,-0.0 0 4 7,1.0 0 1。对 上 述 迭 代 计 算 结 果 进 行 验 证，在 M ATLAB命 令 窗 口 中 输 入 如 下 程 序：A*x输 出 结 果 为：ans=1.0 0 0 00.0 0 0 01.0 0 0 0经 过 验 证，也 确 定 了 解 的 正 确 性。一 般 情 况 下，理 查 森 迭 代 法 难 以 收 敛，因 为 要 使/-/的 谱 半 径 小 于 1,通 常 要 求 力 是 严 格 对 角 占 优 的 矩 阵。197精 通 MATLAB科 学 计 算(第 2 版 i-9.3.3 Jacobi 迭 代 法 如 果 系 数 矩 阵”的 主 对 角 元 全 不 为 0,在 上 节 4 的 分 解 中 取 Q=DC=D-A其 中 O 是 由/的 主 对 角 元 素 组 成 的 对 角 阵，则 有 B=I-D A,r=D b,迭 代 公 式 为：x=(I-D Axk+D b这 种 迭 代 方 法 称 为 Jacobi迭 代 法。在 MATLAB中 编 程 实 现 的 Jacobi迭 代 法 函 数 为：jacobio功 能：用 Jacobi迭 代 法 求 线 性 方 程 组 4 r 的 解 调 用 格 式：x,w=jacobi(A,h,xO,eps,varargin)其 中，A 为 线 性 方 程 组 的 系 数 矩 阵；8 为 线 性 方 程 组 中 的 常 数 向 量；xO为 迭 代 初 始 向 量；e p s为 解 的 精 度 控 制(此 参 数 可 选)；varargin为 迭 代 步 数 控 制(此 参 数 可 选)；x 为 线 性 方 程 组 的 解；为 求 出 所 需 精 度 的 解 实 际 的 迭 代 步 数。Jacobi迭 代 法 的 MATLAB程 序 代 码 如 下：function x,n=jacobi(A,b,xO,eps,varargin)为 采 用 Jacobi迭 代 法 求 线 性 方 程 组 Ax=b的 解 电 线 性 方 程 组 的 系 数 矩 阵：A年 线 性 方 程 组 中 的 常 数 向 量：b卷 迭 代 初 始 向 量：xO优 解 的 精 度 控 制：eps带 迭 代 步 数 控 制：varargin为 线 性 方 程 组 的 解：x生 求 出 所 需 精 度 的 解 实 际 的 迭 代 步 数：nif nargin=3eps=1.Oe-6;M=200;elseif nargin第 9 章 线 性 方 程 组 求 解 M=varargin1;endD=diag(diag(A);务 求 A 的 对 角 矩 阵 L=-tril(Az-1);%求 A 的 下 三 角 阵 U=-triu(A,1);彩 求 A 的 上 三 角 矩 阵 B=D(L+U);f=Db;x=B*xO+f;n=l;%迭 代 次 数 多 迭 代 过 程 while norm(x-xO)=epsxO=x;x=B*xO+f;n=n+l;if(n=M)disp(，Warning:迭 代 次 数 太 多，可 能 不 收 敛！,)；return;endend【例 9-9 Jacobi迭 代 法 求 解 线 性 方 程 组 实 例。用 Jacobi迭 代 法 求 解 以 下 线 性 方 程 组，其 中 取 初 始 值 为 1,1,1。0.9889xj-0.0005x2-0.0002x3=1 x0=ones(3Z1);xz n=jacobi(A,b,xO)输 出 的 计 算 结 果 为：X=1.0114-0.00311.0062输 出 的 迭 代 次 数 为：n=4Y 199精 通 MATLAB科 学 计 算（第 2 版 i-可 见，经 过 4 步 迭 代，Jacobi法 求 出 了 方 程 组 的 解 为：XI,X2,X3=1,0114,-0.0031,1.0062。Jacobi迭 代 法 对 任 意 初 始 向 量 M）收 敛 的 充 要 条 件 为 B 的 谱 半 径 小 于 1,即 8 的 所 有 特 征 值 的 绝 对 值 都 小 于 lo200 第 9 章 线 性 方 程 组 求 解 9.3.4 Gauss-Seidel 迭 代 法 仍 然 设 系 数 矩 阵 N 的 主 对 角 元 全 不 为 零，如 果 对 N 作 以 下 分 解：A=(D-L)-U其 中。的 意 义 同 Jacobi迭 代 法，为 下 三 角 矩 阵，。为 上 三 角 矩 阵，它 的 迭 代 公 式 为:xk+i=(D-)-Uxk+(D-L)-b在 MATLAB中 编 程 实 现 的 Gauss-Seidel迭 代 法 函 数 为：gauseidelo功 能：用 Gauss-Seidel迭 代 法 求 线 性 方 程 组 的 解。调 用 格 式：gauseidel(4 4 的,eps,防。其 中，A 为 线 性 方 程 组 的 系 数 矩 阵；b 为 线 性 方 程 组 中 的 常 数 向 量；xO为 迭 代 初 始 向 量；e p s为 解 的 精 度 控 制(此 参 数 可 选)；M 为 迭 代 步 数 控 制(此 参 数 可 选)；x 为 线 性 方 程 组 的 解；为 求 出 所 需 精 度 的 解 实 际 的 迭 代 步 数。Gauss-Seidel的 MATLAB程 序 代 码 如 下：function x,n=gauseidel(A,b,xO,eps,M)会 采 用 Gauss-Seidel迭 代 法 求 线 性 方 程 组 Ax=b的 解 年 线 性 方 程 组 的 系 数 矩 阵：A%线 性 方 程 组 中 的 常 数 向 量：b他 迭 代 初 始 向 量：xO软 解 的 精 度 控 制：eps多 迭 代 步 数 控 制：M会 线 性 方 程 组 的 解：x务 求 出 所 需 精 度 的 解 实 际 的 迭 代 步 数：nif nargin=3eps=1.Oe-6;M=200;elseif nargin=4M=200;elseif nargin3errorreturn;endD=diag(diag(A);%求 A 的 对 角 矩 阵 201精 通 M ATLAB科 学 计 算（第 2 回 L=-tril(A,-l);U=-triu(A,1);G=(D-L)U；f=(D-L)b；x=G*xO+f;n=l;务 求 A 的 下 三 角 矩 阵%求 A 的 上 三 角 矩 阵 它 迭 代 次 数 国 迭 代 过 程 while norm(x-xO)=epsxO=x;x=G*xO+f;n=n+l;if(n=M)disp(Warning:迭 代 次 数 太 多，可 能 不 收 敛！，)；return;endend【例 9-10】GaussSeidel迭 代 法 求 解 线 性 方 程 组 实 例。用 Gauss Seidel迭 代 法 求 解 以 下 线 性 方 程 组，其 中 取 初 始 值 为 1,1,1。1.4449xi+0.7948x2+0.880 卜 3=1,0.6946XI+1.9568x2+0.1730 x3=00.6213为+0.5226x2+1.9797x3=1 o解：用 Gauss-Seidel迭 代 法 求 解，在 MATLAB命 令 窗 口 中 输 入 求 解 程 序：A=1.4449 0.7948 0.8801;0.6946 1.9568 0.1730;0.6213 0.5226 1.9797;b=l 0 1已；x0=zeros(3,1);x,n=gauseidel(A,b,xO)输 出 的 计 算 结 果 为：X=0.5929-0.24430.3835输 出 的 迭 代 次 数 为：n=n.a可 见，经 过 11步 迭 代，Gauss-Seidel迭 代 法 求 出 方 程 组 的 解 为：Xi,x2,x3=0.5929,-0.2443,0.3835。9.3.5 超 松 弛 迭 代 法 202-第 9 章 线 性 方 程 组 求 解 如 果 对 系 数 矩 阵”作 以 下 分 解：A=(D-a)L)-(l-a)D+a)U)其 中。、L、。的 意 义 与 Gauss-Seidel迭 代 法 相 同。是 一 个 事 先 选 好 的 常 数，称 为 松 弛 因 子，当 时 叫 超 松 弛 法，也 叫 SO R迭 代 法；当。1时 叫 低 松 弛 法。其 迭 代 公 式 为：x*+1=(-尸(1-co)D+a)UXk+(o(D-(oL)b关 于 超 松 弛 迭 代 法 的 收 敛 性 有 如 下 结 论：若 系 数 矩 阵/对 称 正 定，当 0。2,SO R迭 代 法 收 敛。在 MATLAB中 编 程 实 现 的 超 松 弛 迭 代 法 函 数 为：SORo功 能：用 超 松 弛 迭 代 法 求 线 性 方 程 组 的 解。调 用 格 式：x,n=SOR(A,b,xO,w,eps,M)其 中，A 为 线 性 方 程 组 的 系 数 矩 阵；6 为 线 性 方 程 组 中 的 常 数 向 量；xO为 迭 代 初 始 向 量；w 为 松 弛 因 子；e p s为 解 的 精 度 控 制(此 参 数 可 选)；M 为 迭 代 步 数 控 制(此 参 数 可 选)；x 为 线 性 方 程 组 的 解；n 为 求 出 所 需 精 度 的 解 实 际 的 迭 代 步 数。超 松 弛 迭 代 法 的 MATLAB程 序 代 码 如 下：function x,n=SOR(A,b,xO,w,eps,M)带 采 用 超 松 弛 迭 代 法 求 线 性 方 程 组 Ax=b的 解 传 线 性 方 程 组 的 系 数 矩 阵：A省 线 性 方 程 组 中 的 常 数 向 量：b与 迭 代 初 始 向 量：xO%松 弛 因 子：w为 解 的 精 度 控 制：eps超 迭 代 步 数 控 制：M%线 性 方 程 组 的 解：x省 求 出 所 需 精 度 的 解 实 际 的 迭 代 步 数：nif nargin=4eps=1.Oe-6;M=200;elseif nargin 4error 203精 通 MATLAB科 学 计 算(第 2 版 i-returnelseif nargin=5M=200;endif(w=2)与 收 敛 条 件 要 求 error;return;endD=diag(diag(A);%求 A 的 对 角 矩 阵 L=-tril(Az-1);考 求 A 的 下 三 角 矩 阵 U=-triu(Azl);%求 A 的 上 三 角 矩 阵 B=inv(D-L*w)*(1-w)*D+w*U);f=w*inv(D-L*w)*b;x=B*xO+f;n=l;恭 迭 代 次 数 带 迭 代 过 程 while norm(x-xO)=epsx0=x;x=B*xO+f;n=n+l;if(n=M)disp(1 Warning:迭 代 次 数 太 多，可 能 不 收 敛！，)；return;endend【例 9-11 超 松 弛 迭 代 法 求 解 线 性 方 程 组 实 例。用 超 松 弛 迭 代 法 求 解 以 下 线 性 方 程 组，其 中 取 初 始 值 为 0,0,0。4xi+3x2=24 xz n=SOR(Azb,xOz1.25)输 出 的 计 算 结 果 为：X=3.00004.0000-5.0000输 出 的 迭 代 次 数 为：204-第 9 章 线 性 方 程 组 求 解 n=14可 见，经 过 14步 迭 代，SO R迭 代 法 求 出 了 方 程 组 的 解 为 R,X2,X3=3,4,-5。超 松 弛 迭 代 法 还 有 一 种 改 进 形 式，叫 做 对 称 逐 次 超 松 弛 迭 代 法(SSO R),它 采 用 的 是 两 步 迭 代 公 式：(D-a)L)x*+i/2=0(fJxk+b)+(1-a)Dxi；(D-coU)Xk+i=(o(Lxk+b)+(1-v)Z)x*.+i/2在 MATLAB中 编 程 实 现 的 对 称 逐 次 超 松 弛 迭 代 法 函 数 为：SSORo功 能：对 称 逐 次 超 松 弛 迭 代 法 求 线 性 方 程 组 的 解。调 用 格 式：x,=SSOR(4 b,xO,iv,eps,M)o其 中，A 为 线 性 方 程 组 的 系 数 矩 阵；b 为 线 性 方 程 组 中 的 常 数 向 量；xO为 迭 代 初 始 向 量；w 为 松 弛 因 子；e p s为 解 的 精 度 控 制(此 参 数 可 选)；M 为 迭 代 步 数 控 制(此 参 数 可 选)；x 为 线 性 方 程 组 的 解；为 求 出 所 需 精 度 的 解 实 际 的 迭 代 步 数。SSO R算 法 用 MATLAB实 现 如 下：function xz n=SSOR(A,b,xO,w,eps,M)务 采 用 对 称 逐 次 超 松 弛 迭 代 法 求 线 性 方 程 组 Ax=b的 解 金 线 性 方 程 组 的 系 数