版权声明： 本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！

仓库源文，站点原文

title: 数值分析实验 - 线性方程组的直接解法 categories:

数值分析实验 5 - 线性方程组的直接解法

<!-- more -->

目的和意义

给出下列几个不同类型的线性方程组, 请用适当算法计算其解

  1. 线性方程组

    $$ \begin{bmatrix}

      4&2&-3&-1&2&1&0&0&0&0\\
  8&6&-5&-3&6&5&0&1&0&0\\
  4&2&-2&-1&3&2&-1&0&3&1\\
  0&-2&1&5&-1&3&-1&1&9&4\\
  -4&2&6&-1&6&7&-3&3&2&3\\
  8&6&-8&5&7&17&2&6&-3&5\\
  0&2&-1&3&-4&2&5&3&0&1\\
  16&10&-11&-9&17&34&2&-1&2&2\\
  4&6&2&-7&13&9&2&0&12&4\\
  0&0&-1&8&-3&-24&-8&6&3&-1\\
\end{bmatrix}\begin{bmatrix}
  x_1\\x_2\\x_3\\x_4\\x_5\\x_6\\x_7\\x_8\\x_9\\x_{10}
\end{bmatrix}=\begin{bmatrix}
  5\\12\\3\\2\\3\\46\\13\\38\\19\\-21

    \end{bmatrix} $$

  2. 对称正定阵系数阵线方程组

    $$ \begin{bmatrix}

      4&2&-4&0&2&4&0&0\\
  2&2&-1&-2&1&3&2&0\\
  -4&-1&14&1&-8&-3&5&6\\
  0&-2&1&6&-1&-4&-3&3\\
  2&1&-8&-1&22&4&-10&-3\\
  4&3&-3&-4&4&11&1&-4\\
  0&2&5&-3&-10&1&14&2\\
  0&0&6&3&-3&-4&2&19\\
\end{bmatrix}\begin{bmatrix}
  x_1\\x_2\\x_3\\x_4\\x_5\\x_6\\x_7\\x_8
\end{bmatrix}=\begin{bmatrix}
  0\\-6\\20\\23\\9\\-22\\-15\\45

    \end{bmatrix} $$

  3. 三对角形线性方程组

    $$ \begin{bmatrix}

      4&-1&&&&&&&&\\
  -1&4&-1&&&&&&&\\
  &-1&4&-1&&&&&&\\
  &&-1&4&-1&&&&&\\
  &&&-1&4&-1&&&&\\
  &&&&-1&4&-1&&&\\
  &&&&&-1&4&-1&&\\
  &&&&&&-1&4&-1&\\
  &&&&&&&-1&4&-1\\
  &&&&&&&&-1&4\\
\end{bmatrix}\begin{bmatrix}
  x_1\\x_2\\x_3\\x_4\\x_5\\x_6\\x_7\\x_8\\x_9\\x_{10}
\end{bmatrix}=\begin{bmatrix}
  7\\5\\-13\\2\\6\\-12\\14\\-4\\5\\-5

    \end{bmatrix} $$

要求

  1. 对上述三个方程组分别利用 Gauss 顺序消去法与 Gauss 列主元消去法; 平方根法与改进平方根法; 追赶法求解 (选择其一)
  2. 应用结构程序设计编出通用程序
  3. 比较计算结果, 分析数值解误差的原因
  4. 尽可能利用相应模块输出系数矩阵的三角分解式

目的和意义

  1. 通过该课题的实验, 体会模块化结构程序设计方法的优点
  2. 运用所学的计算方法, 解决各类线性方程组的直接算法
  3. 提高分析和解决问题的能力, 做到学以致用
  4. 通过三对角形线性方程组的解法, 体会稀疏线性方程组解法的特点

计算公式

Gauss 消元法

任意情况均可使用

过程即将系数矩阵通过初等行变换化为上三角矩阵或下三角矩阵, 之后回代即可

平方根法

当系数矩阵为对称正定矩阵时可以使用

设方程组为 $Ax=b$

对系数矩阵 $A$ 做分解

$$ A=LL^T $$

其中 $L=(l_{ij})$ 为下三角矩阵,

$$ l{ij}=\begin{cases} \sqrt{a{ii}-\sum{k=1}^{i-1}l{ik}^2},&i=j\ \displaystyle{a{ij}-\sum{k=1}^{j-1}l{ik}l{jk}\over l_{jj}},&i\ne j \end{cases} $$

接下来计算

$$ \begin{cases} Ly=b\ L^Tx=y \end{cases} $$

即可

改进平方根法

当系数矩阵为对称正定矩阵时可以使用

设方程组为 $Ax=b$

对系数矩阵 $A$ 做分解

$$ A=LDL^T $$

其中 $L=(l{ij})$ 为下三角矩阵, $D=(d{i})$ 为对角矩阵, 令 $t{ij}=l{ij}d_j$, 有

$$ t{ij}=a{ij}-\sum{k=1}^{j-1}t{ik}l_{jk} $$

$$ d{i}=a{ii}-\sum{k=1}^{i-1}t{ik}l_{ik} $$

接下来计算

$$ \begin{cases} Ly=b\ L^Tx=D^{-1}y \end{cases} $$

即可

追赶法

当系数矩阵为三对角矩阵时可以使用

设方程组为 $Ax=d$

对系数矩阵做分解

$$ A=\begin{bmatrix} b_1&c_1&&&\ a_2&b_2&c2&&\ &\ddots&\ddots&\ddots&\ &&a{n-1}&b{n-1}&c{n-1}\ &&&a_n&b_n \end{bmatrix}=\begin{bmatrix} \alpha_1&&&&\ \gamma_2&\alpha2&&&\ &\ddots&\ddots&&\ &&\gamma{n-1}&\alpha_{n-1}&\ &&&\gamma_n&\alpha_n \end{bmatrix}\begin{bmatrix} 1&\beta_1&&&\ &1&\beta2&&\ &&\ddots&\ddots&\ &&&1&\beta{n-1}\ &&&&1 \end{bmatrix} $$

有

$$ \alpha_i=b_i-ai\beta{i-1} $$

$$ \beta_i=\frac{c_i}{\alpha_i} $$

$$ \gamma_i=a_i $$

令

$$ L=\begin{bmatrix} \alpha_1&&&&\ \gamma_2&\alpha2&&&\ &\ddots&\ddots&&\ &&\gamma{n-1}&\alpha_{n-1}&\ &&&\gamma_n&\alpha_n \end{bmatrix},~U=\begin{bmatrix} 1&\beta_1&&&\ &1&\beta2&&\ &&\ddots&\ddots&\ &&&1&\beta{n-1}\ &&&&1 \end{bmatrix} $$

接下来计算

$$ \begin{cases} Ly=d\ Ux=y \end{cases} $$

即可

程序设计

主程序

<details open> <summary><font color='orange'>Show code</font></summary> {% icodeweb blog lang:matlab numanaexp-05/main.m %} </details>

输入数据检查

<details open> <summary><font color='orange'>Show code</font></summary> {% icodeweb blog lang:matlab numanaexp-05/input_check.m %} </details>

Gauss 顺序消去法

<details open> <summary><font color='orange'>Show code</font></summary> {% icodeweb blog lang:matlab numanaexp-05/gauss_basic.m %} </details>

Gauss 列主元消去法

<details open> <summary><font color='orange'>Show code</font></summary> {% icodeweb blog lang:matlab numanaexp-05/gauss_advanced.m %} </details>

平方根法

<details open> <summary><font color='orange'>Show code</font></summary> {% icodeweb blog lang:matlab numanaexp-05/sqrt_basic.m %} </details>

改进平方根法

<details open> <summary><font color='orange'>Show code</font></summary> {% icodeweb blog lang:matlab numanaexp-05/sqrt_advanced.m %} </details>

追赶法

<details open> <summary><font color='orange'>Show code</font></summary> {% icodeweb blog lang:matlab numanaexp-05/chase.m %} </details>

结果讨论和分析

结果

分析