发布于  更新于 

数值分析作业 - 线性方程组的直接求解

数值分析

Problem 1

条件数

误差放大因子

Problem 2

Gauss 消元

列主消元

LU 分解

由此前的 Gauss 消元过程

LUP 分解

由列主消元过程

Problem 3

Gauss 消元

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
gaussElimination[A_?SquareMatrixQ, B_?VectorQ] :=
  Module[{a, b, n, x, sum, xmult}, (
    a = A;
    b = B;
    n = Length[A];
    Do[(
      xmult = a[[i, k]]/a[[k, k]];
      a[[i, k]] = xmult;
      Do[(
        a[[i, j]] -= xmult a[[k, j]];
        ), {j, k + 1, n}];
      b[[i]] -= xmult b[[k]];
      ), {k, 1, n - 1}, {i, k + 1, n}];
    x = ConstantArray[0, n];
    x[[n]] = b[[n]]/a[[n, n]];
    Do[(
      sum = b[[i]];
      Do[(
        sum -= a[[i, j]] x[[j]];
        ), {j, i + 1, n}];
      x[[i]] = sum/a[[i, i]];
      ), {i, n - 1, 1, -1}];
    Return[x];
    )];

LU 分解

分解出 LU:

1
2
3
4
5
6
7
8
9
10
11
12
13
luDecomposition[A_?SquareMatrixQ] :=
  Module[{a, n, xmult}, (
    a = A;
    n = Length[A];
    Do[(
      xmult = a[[i, k]]/a[[k, k]];
      a[[i, k]] = xmult;
      Do[(
        a[[i, j]] -= xmult a[[k, j]];
        ), {j, k + 1, n}];
      ), {k, 1, n - 1}, {i, k + 1, n}];
    Return[a];
    )];

利用分解的 LU 和 b 求解 X:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
luSolve[A_?SquareMatrixQ, B_?VectorQ] :=
  Module[{a, b, n, x, sum}, (
    a = A;
    b = B;
    n = Length[A];
    Do[(
      b[[i]] -= a[[i, k]] b[[k]];
      ), {k, 1, n - 1}, {i, k + 1, n}];
    x = ConstantArray[0, n];
    x[[n]] = b[[n]]/a[[n, n]];
    Do[(
      sum = b[[i]];
      Do[(
        sum -= a[[i, j]] x[[j]];
        ), {j, i + 1, n}];
      x[[i]] = sum/a[[i, i]];
      ), {i, n - 1, 1, -1}];
    Return[x];
    )];

分析

1
2
3
4
5
6
a = ( {
    {3, 1, 2},
    {6, 3, 4},
    {3, 2, 5}
   } );
b = {11, 24, 22};

验证 LU 分解的结果:

1
luDecomposition[a] // TeXForm

1
luSolve[luDecomposition[a], b]
1
{1, 2, 3}

验证高斯消元的结果:

1
gaussElimination[a, b]
1
{1, 2, 3}

符号求解 12 阶 Hilbert 矩阵

1
2
3
4
h = HilbertMatrix[12];
b = h . ConstantArray[1, 12];
luSolve[luDecomposition[h], b]
gaussElimination[h, b]
1
2
{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}
{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}

数值求解 12 阶 Hilbert 矩阵:

1
2
3
4
h = HilbertMatrix[12, WorkingPrecision -> MachinePrecision];
b = h . ConstantArray[1, 12];
gaussElimination[a, b] // InputForm
luSolve[luDecomposition[h], b] // InputForm
1
2
3
4
{0.9999999684484674, 1.0000039869936317, 0.9998748290952993,
 1.0017037342786466, 0.9875174669874949, 1.0548274080293765,
 0.847263047196586, 1.2764525211919366, 0.675899878028198,
 1.2373689708614504, 0.9013048730327599, 1.0177833372228902}

数值求解 20 阶 Hilbert 矩阵:

1
2
3
4
h = HilbertMatrix[20, WorkingPrecision -> MachinePrecision];
b = h . ConstantArray[1, 20];
gaussElimination[h, b] // InputForm
luSolve[luDecomposition[h], b] // InputForm
1
2
3
4
5
6
7
{1.0000003421959727, 0.9999417184300158, 1.0024368580495062,
 0.9562365916942744, 1.4191051439028362, -1.3773308447433001,
 9.396245254529573, -17.68936344786112, 26.745626004265336,
 -21.158082774623054, 20.368520692739086, -32.72056597175268,
 47.971136818868764, -25.84867109603557, -15.141625406272288,
 47.602402521150054, -50.83989911260756, 38.52515460862537,
 -15.334245951012626, 4.122978309329297}
许可协议

本文采用 署名-相同方式共享 4.0 国际 许可协议，转载请注明出处。

分享文章
快来参与讨论吧