发布于  更新于 

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

数值分析

Problem 1

20230324202737
20230324202737

条件数

误差放大因子

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}