线性方程组直接解法

一开始还是有点虚的，NumPy 和 Pandas 都是前学后忘，于是第一次作业拖到周末才写，意外地顺利。

反正就是不停地写循环，略烦。

\(LDL^T\) 法

算法就是 \(LDL^T\) 法。书上啥都有。

比较痛的一点是，这个方法为了提高效率把 \(L\) 和 \(D\) 直接存储在 \(A\) 里了，但是最后计算误差的时候需要调用原值 \(A\)，所以做分解之前应该先备份一个原矩阵。向量 \(b\) 的备份同理。

实现如下：

# 任务一：改进的平方根法求解线性方程组

import numpy as np
from numpy.linalg import norm

n = 32
H = 1./(np.arange(1, n + 1)+np.arange(0, n)[:,np.newaxis])  #构造n阶Hilbert矩阵
xtrue = np.ones((n, 1))
x = np.ones((n, 1))
b = H.dot(xtrue)

#待定系数法计算LDL^t分解
K = H.copy()
c = b.copy()
for j in range(0, n):
    for k in range(0, j):
        K[j, j] = K[j, j] - K[j, k] * K[j, k] * K[k, k]
    for i in range(j + 1, n):
        for k in range(0, j):
            K[i, j] = K[i, j] - K[i, k] * K[j, k] * K[k, k]
        K[i,j] = K[i, j] / K[j, j]

#写出L和D的具体形式
L = np.zeros((n, n))
D = np.zeros((n, n))
for i in range(0, n):
    L[i, i] = 1
    for j in range(0, i):
        L[i, j] = K[i, j]
    D[i, i] = K[i, i]

#解方程Ly=b
y = np.zeros((n, 1))
y[0] = c[0]
for i in range(1, n):
    for j in range(0, i):
        c[i] = c[i] - K[i, j] * y[j]
    y[i] = c[i]

#解方程DL^tx=y
x = np.zeros((n, 1))
x[n-1] = c[n-1] / K[n-1,n-1]
for i in range(n-2, -1, -1):
    x[i] = y[i] / K[i,i]
    for j in range(i+1, n):
        x[i] = x[i] - K[j,i] * x[j]

#输出结果
print('相对残量为: %.4f'%(norm(b - np.matmul(H, x)) / norm(b)))
print('相对误差为: %.4f'%(norm(x - xtrue) / norm(xtrue)))

输出结果：

1 2	相对残量为: 0.0000 相对误差为: 111.3149

残量为 \(0\) 说明分解正确。考虑到 \(Hilbert\) 矩阵是个病态矩阵，这个值还算可以接受。虽然真值是一个分量全为 \(1\) 的 \(32\) 阶向量，而数值解从第 \(7\) 个分量开始就大量出现绝对值是两位数的情况了。

\(LDL^T\) 分解从结果来看仍然是 \(LU\) 分解的一种变形，但在本实验中，利用 \(LU\) 分解得到的数值解各分量量级偏小，甚至达到 \(10^{-5}\) 级别。输出结果为：

1 2	相对残量为: 0.0000 相对误差为: 73385.0585

相对误差大了两个量级，说明对于对称矩阵来说 \(LDL^T\) 分解是更好的选择。

我在做数值实验的过程中先解决了 \(LU\) 分解，回头又做了 \(LDL^T\) 分解，因此二者的结果不能很好地匹配也造成了不小的困难..

\(LU\) 法

我感觉最难的是从 .mat 文件里载入矩阵...

\(LU\) 法，学过线代就能做，我大一上学期的 C++ 大作业里还写这个了..
实现如下：

# 任务二：实现按列存储的JKI型LU分解

from scipy.io import loadmat
import numpy as np
from numpy.linalg import norm

def LU_JKI(B):
    A = B.copy()
    (m, n) = A.shape  #本任务中m=n
    L = np.zeros((m, n))
    U = np.zeros((m, n))
    #计算L矩阵和U矩阵

    for i in range(1, n):
        for k in range(0, i):
            A[i, k] = A[i, k] / A[k, k]
            for j in range(k + 1, n):
                A[i, j] = A[i, j] - A[i, k] * A[k, j]
    
    for i in range(0, n):
        L[i, i] = 1
        for j in range(0, i):
            L[i, j] = A[i, j]
        for j in range(i, n):
            U[i, j] = A[i, j]

    return L,U

#主函数部分
A = loadmat('./data.mat')['A'] #载入数据

L,U = LU_JKI(A)

#测试，计算相对误差 ||A - L*U|| / ||A||，这里用 F-范数
print('relerr = %.4f'%(norm(A - np.matmul(L, U),'fro') / norm(A, 'fro')))

类似地，注意算法虽然是在 \(A\) 上进行的，但在程序中需要保留 \(A\) 的原值来计算相对误差，因此在函数中使用的应该是 \(A\) 的复制，也就有了 A=B.copy() 这一步。

输出结果：

1	relerr = 0.0000

相对误差接近于 \(0\)，说明实验是成功的。

后记

不是，这玩意不是我后来咕咕了没再写，而是布置过一次作业之后就再也没有了（

本来说是基本每一讲后面都有数值实验的，结果不知道是老师咕咕还是助教咕咕，总之是并没有再布置过实验，少了很多乐趣。之前看过贵系数值分析的资料 Repo，两边的风格还是很不一样的。

就这么草率地结束了..

『姑妄言之姑妄听之』

数值实验顺便 Python 复健

线性方程组直接解法

\(LDL^T\) 法

\(LU\) 法

后记