ComPhys-2024 – Chapt 8: Numerical methods for matrices

Show the code

import numpy as np
import matplotlib.pyplot as plt
from sympy import *
from IPython.display import Markdown
from sympy.solvers.ode.systems import dsolve_system
from sympy import *
import array_to_latex as a2l

很多数值以及解析物理问题中都涉及到矩阵方程，因此矩阵方程的求解方法可以应用到相关的ODE以及PDE中。例如，PDE形式的EVP可以重写程矩阵形式，BVP问题离散后也可以写成一系列的线性代数方程。

1 物理中的矩阵问题

1.1 分子振动能级

物理中很多问题都可以写成矩阵形式，例如，我们想研究一个具有 $n$ 个振动自由度分子的振动谱，首先就是在平衡结构附近将系统势能展开成广义坐标的二阶导数，进而研究系统的简谐振动：

\[ U(q_1,q_2,\dots,q_n) \simeq \frac{1}{2}\sum_{i,j=1}^nA_{ij}q_iq_j \]

其中 $q_i$ 时广义坐标， $A_{ij}$ 时广义弹性常数矩阵元（可以通过量子化学计算或者实验获得）。这里将平衡势能取为 $0$ 点，系统的动能一般可以写为

\[ T(\dot{q}_1, \dot{q}_2, \dots, \dot{q}_n) \simeq \frac{1}{2}\sum_{i,j=1}^nM_{ij}\dot{q}_i\dot{q}_j, \]

其中 $\dot{q}_i = dq_i/dt$ 为广义速度， $M_{ij}$ 为广义质量矩阵元，其值取决于分子种类。应用拉格朗日方程：

\[ \frac{\partial \mathcal{L}}{\partial t} - \frac{d}{dt} \frac{\partial \mathcal{L}}{\partial \dot{q}_i} = 0, \]

其中 $\mathcal{L} = T-U$ 为系统的拉格朗日量：

\[ \sum_{j=1}^n(A_{ij}q_j+M_{ij}\ddot{q}_j) = 0, \]

其中 $i = 1, 2, \dots, n$. 假设广义坐标振荡性地依赖于时间：

\[ q_j = x_j e^{-i\omega t}, \]

代入后有

\[ \sum_{j = 1}^n(A_{ij} - \omega^2M_{ij})x_j = 0, \]

可以写成矩阵形式：

\[ \begin{pmatrix} A_{11} & \dots & A_{1n} \\ \vdots & \vdots & \vdots \\ A_{n1} & \dots & A_{nn} \end{pmatrix} \begin{pmatrix} x_1 \\ \vdots \\ x_n \end{pmatrix} = \lambda \begin{pmatrix} M_{11} & \dots & M_{1n} \\ \vdots & \vdots & \vdots \\ M_{n1} & \dots & M_{nn} \end{pmatrix} \begin{pmatrix} x_1 \\ \vdots \\ x_n \end{pmatrix} \]

或者

\[ \mathbf{A x} = \lambda \mathbf{M x}, \]

其中 $\lambda = \omega^2$ 为本征值，$\mathbf{x}$ 为相应的本征矢量。此方程为一系列的齐次方程组，为了得到非零解，系数矩阵的行列式必须为0 \[ |\mathbf{A} - \lambda \mathbf{M}| = 0 \] 该久期方程的根$\lambda_k, k = 1, 2, \dots, n$ 给出了原子所有可能的振动频率 $\omega_k = \sqrt{\lambda_k}$。

1.2 电路

通过基尔霍夫定律，我们可以得到一系列的电压、电流方程，进而去求解其他未知量。我们以惠斯通电桥(Wheatstone bridge)为例，其中有三个独立的回路，我们选择第一个为穿过源和两个上桥；第二和第三个为电流表的左边和右边的回路。每个回路都可以得到一个方程： \[ \begin{aligned}[r] r_s i_1 + r_1 i_2 + r_2 i_3 = v_0, \\ -r_x i_1 + (r_1 + r_x + r_a)i_2 - r_ai_3 = 0, \\ -r_3 i_1 - r_a i_2 + (r_2 + r_3 + r_a)i_3 = 0, \end{aligned} \]

上述方程组可以写成矩阵形式 \[ \mathbf{R i} = \mathbf{v}, \] 其中 \[ \mathbf{R} = \begin{pmatrix} r_s & r_1 & r_2 \\ -r_x & r_1 + r_x + r_a & -r_a \\ -r_3 & -r_a & r_2 + r_3 + r_a \end{pmatrix} \] 为阻抗系数矩阵， \[ \mathbf{i} = \begin{pmatrix} i_1 \\ i_2 \\ i_3 \end{pmatrix} ~~~~ \text{and} ~~~~ ~~~~ \mathbf{v} = \begin{pmatrix} v_0 \\ 0 \\ 0 \end{pmatrix} \] 为电流和电压矩阵，方程两端都乘以逆矩阵$\mathbf{R}^{-1}$ \[ \mathbf{i} = \mathbf{R}^{-1} \mathbf{v}, \]

1.3 多电子系统电子结构

这里以 $H_3^+$ 为例，我们可以得到一个非常有用的系统。三个质子组成了一个等边三角形，并且共享了系统中的两个电子。假设我们可以用一个简单的哈密顿量 $\mathcal{H}$ 来描述这个系统，包含一项描述电子从一个位置跃迁到另一个位置，以及另一项用于描述双占据位置上两个电子之间的排斥： \[ \mathcal{H} = -t\sum_{i \neq j} a^\dagger_{i \sigma} a_{j \sigma} + U \sum n_{i\uparrow} n_{i\downarrow}, \] 其中$a^\dagger_{i \sigma}$和$a_{j \sigma}$代表自旋向上$(\sigma = \uparrow)$ 和向下$(\sigma = \downarrow)$的产生和湮灭算符。 $n_{i\sigma} = a^\dagger_{i \sigma}a_{j \sigma}$是第 $i$ 层上的占据数。当用于描述高度关联电子体系时，这个哈密顿量称为Hubbard模型。其中的$t$和$U$可以从量子化学或者实验中得到。系统的Schro"dinger方程为 \[ \mathcal{H}\ket{\Psi_k} = \mathcal{E}_k \ket{\Psi_k}. \] 因为每个位置只有一个轨道，所以在单粒子表示下，两个电子总共有15个可能的状态。 \[ \begin{aligned}[l] \ket{\phi_1} = a^\dagger_{1\uparrow} a^\dagger_{1\downarrow} \ket{0}, \ket{\phi_2} = a^\dagger_{2\uparrow} a^\dagger_{2\downarrow} \ket{0}, \ket{\phi_3} = a^\dagger_{3\uparrow} a^\dagger_{3\downarrow} \ket{0}, \\ \ket{\phi_4} = a^\dagger_{1\uparrow} a^\dagger_{2\downarrow} \ket{0}, \ket{\phi_5} = a^\dagger_{1\uparrow} a^\dagger_{3\downarrow} \ket{0}, \ket{\phi_6} = a^\dagger_{2\uparrow} a^\dagger_{3\downarrow} \ket{0}, \\ \ket{\phi_7} = a^\dagger_{1\downarrow} a^\dagger_{2\uparrow} \ket{0}, \ket{\phi_8} = a^\dagger_{1\downarrow} a^\dagger_{3\uparrow} \ket{0}, \ket{\phi_9} = a^\dagger_{2\downarrow} a^\dagger_{3\uparrow} \ket{0}, \\ \ket{\phi_10} = a^\dagger_{1\uparrow} a^\dagger_{2\uparrow} \ket{0}, \ket{\phi_11} = a^\dagger_{1\uparrow} a^\dagger_{3\uparrow} \ket{0}, \ket{\phi_12} = a^\dagger_{2\uparrow} a^\dagger_{3\uparrow} \ket{0}, \\ \ket{\phi_13} = a^\dagger_{1\downarrow} a^\dagger_{2\downarrow} \ket{0}, \ket{\phi_14} = a^\dagger_{1\downarrow} a^\dagger_{3\downarrow} \ket{0}, \ket{\phi_15} = a^\dagger_{2\downarrow} a^\dagger_{3\downarrow} \ket{0}, \end{aligned} \] 其中$\ket{0}$为真空态。此时，Schrodinger方程可以写成矩阵形式： \[ \mathbf{H \Psi}_k = \mathcal{E}_k \mathbf{\Psi}_k \] 其中$H_{ij} = \bra{\phi_i}\mathcal{H}\ket{\phi_j}$，其中$i, j = 1, 2, \dots, 15$，$\Psi_{ki} = \braket{\phi_i|\Psi_k}$，最后也可以归结为一个久期方程： \[ |\mathbf{H} - \mathcal{E}I| = 0 \] 其中$\mathbf{I}$ 为单位矩阵 $I_{ij} = \delta_{ij}$。这里的量子问题与分子的经典振动具有相似的数学形式。我们可以利用系统对称性对问题进行简化。总自旋以及总自旋的$z$分量与哈密顿量对易，因此，它们是好量子数。

此时，我们可以将哈密顿矩阵越化为最小块大小为$2 \times 2$的块对角矩阵。

当得到系统所有的本征值和本征矢量，就可以分析$H_3^+$的电子、光学、磁学性质。

2 Basic matrix operations

\[ \mathbf{Ax} = \mathbf{b} \leftrightarrow \sum_{j=1}^nA_{ij}x_j = b_i, \]
\[ \mathbf{C} = \mathbf{AB} \leftrightarrow C_{ij} = \sum_{k=1}^nA_{ik}B_{kj}, \]
\[ \mathbf{AA}^{-1} = \mathbf{A}^{-1}\mathbf{A} = \mathbf{I} \]
Derminant: \[ |\mathbf{A}| = \sum_{i=1}^n (-1)^{i+j}A_{ij}|\mathbf{R}_{ij}| \]

$|\mathbf{R}_{ij}|$ is the determinant of the residual matrix $\mathbf{R}_{ij}$ of $\mathbf{A}$ with its $i$th row and $j$th column removed.

$C_{ij} = (-1)^{i+j}|\mathbf{R}_{ij}|$ is called a cofactor of $\mathbf{A}$.

Inverse: \[ A^{-1}_{ij} = \frac{C_{ij}}{|\mathbf{A}|}, \]
If $|\mathbf{A}| \neq 0$, it is called a nonsingular matrix.
If $|\mathbf{A}| = 0$, it is called a singular matrix.
Trace of $\mathbf{A}$: sum of all its diagonal elements, $tr(\mathbf{A}) = \sum_{i=1}^n A_{ii}$.
Transpose: \[ A^T_{ij} = A_{ji} \]
Orthogonality: \[ \mathbf{A}^T\mathbf{A} = \mathbf{AA}^T = \mathbf{I} \]
Hermitian: complex conjugate $\mathbf{A}^\dagger \equiv \mathbf{A}^{T*}$ \[ \mathbf{A}^\dagger = \mathbf{A}^{-1} \]
Add one row/column multiplied by a factor $\lambda$ to another row/column \[ A'_{ij} = A_{ij} + \lambda A_{kj}, \] This operation preserves the determinant, $|\mathbf{A}'| = |\mathbf{A}|$ and can be represented by a matrix multiplication \[ \mathbf{A}' = \mathbf{M}\mathbf{A} \] 其中，$\mathbf{M}$具有单位对角元，再加上一个非0非对角元 $M_{ik} = \lambda$。当调换任意两个行或者列，行列式值会变号。

$H^+_3$电子结构可以定义为 \[ \mathbf{Ax} = \lambda \mathbf{x} \] 其中 $\mathbf{x}$ 和 $\lambda$ 为矩阵的本征矢量以及相应的本征值。

如果将 $\mathbf{M}^{-1}\mathbf{A} = \mathbf{B}$ 作为问题的矩阵，$\mathbf{Ax} = \lambda \mathbf{Mx}$，将有 $\mathbf{Bx} = \lambda \mathbf{x}$。

矩阵本征值问题可以看作是一个迭代求解的线性方程组问题。例如如果我们想求解上述方程的本征值和本征矢量，可以使用递归法 \[ \mathbf{Ax}^{(k+1)} = \lambda^{(k)} + \mathbf{x}^{(k)} \] 其中$\lambda^{(k)}$和$\mathbf{x}^{(k)}$ 为第$k$次迭代中合适的本征值以及本征矢量，即，迭代法。

当我们选择$\mathbf{x} = \mathbf{Sy}$时，由于 \[ \mathbf{By} = \lambda\mathbf{y} \] 与 \[ \mathbf{By} = \lambda\mathbf{y} \] 等价，矩阵本征值在相似变换 \[ \mathbf{B} = \mathbf{S}^{-1}\mathbf{A}\mathbf{S} \] 下与矩阵$\mathbf{A}$本征值之一相同，

其中$\mathbf{S}$矩阵非奇异，同时也意味着 \[ |\mathbf{B}| = |\mathbf{A}| = \prod_{i=1}^n \lambda_i \]

3 线性方程系统

对角线以上、以下矩阵元为0的矩阵成为上三角、下三角矩阵。求解矩阵最简单的方法就是高斯消元法，就是首先使用基本矩阵操作将系数矩阵变换成上、下三角矩阵，再通过向后、前替换求解方程的解。矩阵的逆以及行列式也可以通过这种方式求的。

以\[ \mathbf{Ax} = \mathbf{b}\] 为例。如果$|\mathbf{A}| \neq 0$以及$\mathbf{b} \neq 0$，系统就具有唯一解。将原来的矩阵标记为$\mathbf{A}^(0)$，$\mathbf{A}^(j)$为$j$次操作后的矩阵。$\mathbf{b}$也采取同样的标记。

3.1 高斯消元

高斯消元的基本思想是将原来的矩阵变换为具有相同解的上、下三角矩阵，这里以变换成上三角为例。同样也可以通过交换行与列的角色来得到下三角矩阵。在每次操作中，我们消掉不属于上三角的列中元素。最终的线性方程组为 \[ \mathbf{A}^{(n-1)}\mathbf{x} = \mathbf{b}^{(n-1)} \] 其中$A_{ij}^{(n-1)} = 0~~\text{for}~~i>j$，这个过程很简单：

Gauss elimination

首先给第一个方程乘以$-A_{i1}^{(0)}/A_{11}^{(0)}$，也就是矩阵的第一行以及$b_1^{(0)}$，并将其加到第$i$个方程对于所有的$i>1$；
除第一行外，矩阵中每行的第一个元素都被消去，矩阵此时变成$\mathbf{A}^{(1)}$；
第二个方程乘以$-A_{i2}^{(1)}/A_{22}^{(1)}$，并将其加到其他方程$i>2$；
除第一、二行的第二个元素，其他行的第二个元素都别消掉了；
对第三、四、$n-1$个方程重复该过程，该系数矩阵就变成了上三角矩阵 $\mathbf{A}^{(n-1)}$；

上、下三角系数矩阵可以通过向后代入轻易求解。

由于所有的对角线元都会被用到分母中，当任意一个分母为0时，该方法将不再适用。这可以通过行调换解决：将 $k$ 行与 $k$ 行以下 $k$ 列中最大元素 $A_{rk}$ 所在行$r$行进行调换。该过程只是调换了各个方程的顺序，不改变方程的解。

这里我们考虑部分调换情况，即只对特定列的剩余行进行调换。完全调换方法中会在所有剩余元素中搜索调换元。这里以部分调换为例。

我们先从$|A^{(0)}_{1i}|$中搜索最大的元素，假设为$A^{(0)}_{k_1 1}$，之后将第一行与$k_1$行调换，从而消除第一行以外每行的第一个元素。之后对每行都进行类似操作就可以将$A$矩阵变换成上三角矩阵。

物理上，在寻找枢轴元时我们不需要交换矩阵中的行。可以用一个指标来记录所有行中枢轴元的顺序。此外，对于$i>j$时，我们需要$A^{(j-1)}_{ij} / A^{(j-1)}_{k_j j}$比率来调整$\mathbf{b}^{(j-1)}$到$\mathbf{b}^j$，并存储到与$\mathbf{A}^{(n-1)}$相同二维矩阵中$A_{ij}^{(n-1)} = 0$的位置，即，对角线以下。当我们确定枢轴元时，为了得到公平的比较时，我们需要将每行的元素用所在行最大的元素进行重新标度。这个重新标度同样减少了潜在的舍入误差。

Show the code

import numpy as np
def gauss_elim(B):
    A = B.copy()
    nn = len(A)
    sig = 1.0
    for i in range(nn):
        if (A[i, i] != 0):
            for j in range(i+1, nn):
                A[j, :] = A[j, :] - A[j, i] * A[i, :] / A[i, i]
        else:
            index = np.where(np.abs(A[i:,i]) == np.max(np.abs(A[i:,i])))
            temp = A[i, :].copy()
            A[i, :] = A[index[0][0]+i, :].copy()
            A[index[0][0]+i, :] = temp
            sig = sig * -1.0
            if (A[i, i] == 0):
                continue
            for j in range(i+1, nn):
                A[j, :] = A[j, :] - A[j, i] * A[i, :] / A[i, i]
    return A, sig

AA = np.array([[1, 2, 3], [4, 8, 6], [7, 8, 10]], dtype=float)
BB, sig = gauss_elim(AA)
print('A\n', AA)
print("A'\n", BB)

A
 [[ 1.  2.  3.]
 [ 4.  8.  6.]
 [ 7.  8. 10.]]
A'
 [[  1.   2.   3.]
 [  0.  -6. -11.]
 [  0.   0.  -6.]]

3.2 矩阵行列式

得到上三角或者下三角后，就很容易得到行列式。高斯消元只改变行列式的符号，不改变值，而符号可以通过行列的轮换次数确定。对于上三角矩阵，其行列式由对角线元素乘积给出。

Show the code

res = 1.0
for i in range(len(BB)):
    res = res * BB[i, i]
print(res * sig)

print(np.linalg.det(AA))

print(np.linalg.det(BB) * sig)

-36.0
-36.0
-36.0

3.3 线性方程组解

类似的，通过高斯消元法将矩阵进行对角化的方法也可以用来求解线性方程组，解可以通过向后代入求的： \[ x_i = \frac{1}{A^{(n-1)}_{k_i i}} \left( b^{(n-1)}_{k_i} - \sum_{j=i+1}^n A^{(n-1)}_{k_i j} x_j \right) \] 对于 $i = n-1, n - 2, \dots, 1$，出发点为 $x_n = b_{k_n}^{(n-1)} / A_{k_n n}^{(n-1)}$。这里使用 $k_i$ 代表 $i$ 列枢轴元素的行指标。容易证实，上述递归满足线性方程组 $\mathbf{A}^{(n-1)}\mathbf{x} = \mathbf{b}^{(n-1)}$。这里我们以惠斯通电桥中求解电流为例，其中 $r_s = r_1 = r_2 = r_x = r_a = 100~\Omega$ 以及 $v_0 = 200~\text{V}$。

Show the code

AA = np.array([[100.0, 100.0, 100.0], [-100.0, 300.0, -100.0], [-100.0, -100.0, 300.0]])
b = np.array([200.0, 0.0, 0.0]).reshape(3, 1)
problem = np.concatenate((AA, b), axis = 1)
print(problem.shape)
print(problem)

res, sig = gauss_elim(problem)

A = res[:, 0:3]
B = res[:, -1]


x = np.zeros(len(AA))
x[-1] = B[-1] / A[-1, -1]
for n in range(len(AA)-2, -1, -1):
    sm = 0.0
    for i in range(n+1, len(AA)):
        sm += A[n, i] * x[i]
    x[n] = 1.0 / A[n, n] * (B[n] - sm)
print(x)

(3, 4)
[[ 100.  100.  100.  200.]
 [-100.  300. -100.    0.]
 [-100. -100.  300.    0.]]
[1.  0.5 0.5]

3.4 矩阵的逆

$\mathbf{A}$ 的逆矩阵 $\mathbf{A}^{-1}$ 可以通过类似方式求的： \[ A^{-1}_{ij} = x_{ij}, \] 其中 $x_{ij}$ 为方程 $\mathbf{A}\mathbf{x}_j = \mathbf{b}$ 的解，且，$b_{ij} = \delta_{ij}$。如果将 $\mathbf{b}$ 展开成单位矩阵来求解线性方程，单位矩阵每列的解组成了 $\mathbf{A}^{-1}$ 的逆矩阵，此时只需要对上面程序进行稍加修改即可。

Show the code

def solveeq(AA, BB):
    problem = np.concatenate((AA, BB), axis = 1)
    res, sig = gauss_elim(problem)

    A = res[:, 0:3]
    B = res[:, -1]

    x = np.zeros(len(A))
    x[-1] = B[-1] / A[-1, -1]
    for n in range(len(A)-2, -1, -1):
        sm = 0.0
        for i in range(n+1, len(A)):
            sm += A[n, i] * x[i]
        x[n] = 1.0 / A[n, n] * (B[n] - sm)
    return x

AA = np.array([[100.0, 100.0, 100.0], [-100.0, 300.0, -100.0], [-100.0, -100.0, 300.0]])

def inverse(AA):
    result = np.zeros(AA.shape)
    for i in range(len(AA)):
        BB = np.zeros(len(AA)).reshape(len(AA), 1)
        BB[i] = 1.0
        result[:, i] = solveeq(AA, BB)
    return result

result = inverse(AA)

print(AA)
print(result)

print(np.dot(AA, result))

print(np.linalg.inv(AA))

[[ 100.  100.  100.]
 [-100.  300. -100.]
 [-100. -100.  300.]]
[[ 0.005  -0.0025 -0.0025]
 [ 0.0025  0.0025  0.    ]
 [ 0.0025  0.      0.0025]]
[[ 1.00000000e+00  5.20417043e-18  5.20417043e-18]
 [-5.20417043e-18  1.00000000e+00 -5.20417043e-18]
 [ 1.56125113e-17 -5.20417043e-18  1.00000000e+00]]
[[ 0.005  -0.0025 -0.0025]
 [ 0.0025  0.0025  0.    ]
 [ 0.0025  0.      0.0025]]

3.5 LU分解

与高斯消元法相关的方法是LU分解，其将矩阵分解成一个上三角和一个下三角矩阵的乘积。如前所属，高斯消元法对应于$n-1$步矩阵操作并且可以表示为矩阵乘法 \[ \mathbf{A}^{(n-1)} = \mathbf{MA}, \qquad(1)\] 其中 \[ \mathbf{M} = \mathbf{M}^{(n-1)}\mathbf{M}^{(n-2)}\dots\mathbf{M}^{(1)}, \] 其中每个 $\mathbf{M}^{(r)},~~\text{for}~~ r = 1, 2, \dots, n-1$，代表一步消元。容易验证 $\mathbf{M}$ 是下三角矩阵其对角元为 $-1$, $M_{k_ij} = -A_{k_i j}^{(j-1)} / A_{k_j j}^{(j-1)}~~~\text{for}~~~ i > j$，方程(Equation 1)等价于 \[ \mathbf{A} = \mathbf{LU}. \] 其中 $\mathbf{U} = \mathbf{A}^{(n-1)}$ 为上三角矩阵，$\mathbf{L} = \mathbf{M}^{-1}$ 为下三角矩阵。注意下三角矩阵的逆矩阵仍然为一个下三角矩阵。由于 $M_{ii} = -1$，$L_{ii} = 1$，$L_{ij} = -M_{ij} = A_{k_ij}^{(j-1)}/A_{k_jj}^{(j-1)} ~~~ \text{for} ~~~ i>j$，这就是之前用于部分轮换高斯消元时存储的比例。这个 $\mathbf{LU}$ 分解的方法将 $\mathbf{U}$ 的非 0 元素以及 $\mathbf{L}$ 的非0，非对角元素存储到了返回矩阵。

LU分解中选择 $L_{ii} = 1$ 的方法称为 Doolittle 分解。另外一个方法是选择 $U_{ii} = 1$；这称为Crout分解。Crout分解等价于用 $U_{ii}$ 将 $U_{ij}$ 缩放对于 $i<j$，以及给 $L_{ij}$ 乘以 $U_{jj}$ 对于 $i>j$（ $L_{ij}$ 来自于Doolitte分解）。

一般情况下，$\mathbf{L}$ 和 $\mathbf{U}$ 的元素可以从二者乘积与 $\mathbf{A}$ 对比得到，结果可以写成一个迭代形式： \[ \begin{aligned}[l] L_{ij} = \frac{1}{U_{ij}}\left(A_{ij} - \sum_{k=1}^{j-1}L_{ik}U_{kj}\right), \\ U_{ij} = \frac{1}{L_{ij}}\left(A_{ij} - \sum_{k=1}^{i-1}L_{ik}U_{kj}\right), \end{aligned} \] 其中 $L_{i1} = A_{i1}/U_{11}, U_{1j} = A_{1j} / L_{11}$ 为起始点。以此作为迭代求解，仍然需要考虑轮换的问题。

--- title: "Chapt 8: Numerical methods for matrices" from: markdown+emoji navbar: false language: en: title-block-author-single: "Author" title-block-published: "Updated" zh: title-block-author-single: "作者" title-block-published: "更新时间" author: - name: "Feng Wan" email: wanfeng@xjtu.edu.cn orcid: 0000-0002-9715-3032 affiliations: - name: "B832@仲英楼 or 3133@泓理楼 // School of Physcis @ Xian Jiaoton University" address: "B832@仲英楼 or 3133@泓理楼" --- ```{=html} <style> red { color: red } yellow { color: yellow } blue {color: blue} </style> ``` ```{python} import numpy as np import matplotlib.pyplot as plt from sympy import * from IPython.display import Markdown from sympy.solvers.ode.systems import dsolve_system from sympy import * import array_to_latex as a2l ``` > 很多数值以及解析物理问题中都涉及到矩阵方程，因此矩阵方程的求解方法可以应用到相关的ODE以及PDE中。例如，PDE形式的EVP可以重写程矩阵形式，BVP问题离散后也可以写成一系列的线性代数方程。 ## 物理中的矩阵问题 ### 分子振动能级物理中很多问题都可以写成矩阵形式，例如，我们想研究一个具有 \$n\$ 个振动自由度分子的振动谱，首先就是在平衡结构附近将系统势能展开成广义坐标的二阶导数，进而研究系统的简谐振动： $$ U(q_1,q_2,\dots,q_n) \simeq \frac{1}{2}\sum_{i,j=1}^nA_{ij}q_iq_j $$ 其中 $q_i$ 时广义坐标， $A_{ij}$ 时广义弹性常数矩阵元（可以通过量子化学计算或者实验获得）。这里将平衡势能取为 $0$ 点，系统的动能一般可以写为 $$ T(\dot{q}_1, \dot{q}_2, \dots, \dot{q}_n) \simeq \frac{1}{2}\sum_{i,j=1}^nM_{ij}\dot{q}_i\dot{q}_j, $$ 其中 $\dot{q}_i = dq_i/dt$ 为广义速度， $M_{ij}$ 为广义质量矩阵元，其值取决于分子种类。应用拉格朗日方程： $$ \frac{\partial \mathcal{L}}{\partial t} - \frac{d}{dt} \frac{\partial \mathcal{L}}{\partial \dot{q}_i} = 0, $$ 其中 $\mathcal{L} = T-U$ 为系统的拉格朗日量： $$ \sum_{j=1}^n(A_{ij}q_j+M_{ij}\ddot{q}_j) = 0, $$ 其中 $i = 1, 2, \dots, n$. 假设广义坐标振荡性地依赖于时间： $$ q_j = x_j e^{-i\omega t}, $$ 代入后有 $$ \sum_{j = 1}^n(A_{ij} - \omega^2M_{ij})x_j = 0, $$ 可以写成矩阵形式： $$ \begin{pmatrix} A_{11} & \dots & A_{1n} \\ \vdots & \vdots & \vdots \\ A_{n1} & \dots & A_{nn} \end{pmatrix} \begin{pmatrix} x_1 \\ \vdots \\ x_n \end{pmatrix} = \lambda \begin{pmatrix} M_{11} & \dots & M_{1n} \\ \vdots & \vdots & \vdots \\ M_{n1} & \dots & M_{nn} \end{pmatrix} \begin{pmatrix} x_1 \\ \vdots \\ x_n \end{pmatrix} $$ 或者 $$ \mathbf{A x} = \lambda \mathbf{M x}, $$ 其中 $\lambda = \omega^2$ 为本征值，$\mathbf{x}$ 为相应的本征矢量。此方程为一系列的齐次方程组，为了得到非零解，系数矩阵的行列式必须为0 $$ |\mathbf{A} - \lambda \mathbf{M}| = 0 $$ 该久期方程的根$\lambda_k, k = 1, 2, \dots, n$ 给出了原子所有可能的振动频率 $\omega_k = \sqrt{\lambda_k}$。 ### 电路 ::::{.columns} :::{.column width="50%"} 通过基尔霍夫定律，我们可以得到一系列的电压、电流方程，进而去求解其他未知量。我们以惠斯通电桥(Wheatstone bridge)为例，其中有三个独立的回路，我们选择第一个为穿过源和两个上桥；第二和第三个为电流表的左边和右边的回路。每个回路都可以得到一个方程： $$ \begin{aligned}[r] r_s i_1 + r_1 i_2 + r_2 i_3 = v_0, \\ -r_x i_1 + (r_1 + r_x + r_a)i_2 - r_ai_3 = 0, \\ -r_3 i_1 - r_a i_2 + (r_2 + r_3 + r_a)i_3 = 0, \end{aligned} $$ ::: :::{.column width="10%"} ::: :::{.column width="40%"} ![](figures/wheatstone.png) ::: :::: --- 上述方程组可以写成矩阵形式 $$ \mathbf{R i} = \mathbf{v}, $$ 其中 $$ \mathbf{R} = \begin{pmatrix} r_s & r_1 & r_2 \\ -r_x & r_1 + r_x + r_a & -r_a \\ -r_3 & -r_a & r_2 + r_3 + r_a \end{pmatrix} $$ 为阻抗系数矩阵， $$ \mathbf{i} = \begin{pmatrix} i_1 \\ i_2 \\ i_3 \end{pmatrix} ~~~~ \text{and} ~~~~ ~~~~ \mathbf{v} = \begin{pmatrix} v_0 \\ 0 \\ 0 \end{pmatrix} $$ 为电流和电压矩阵，方程两端都乘以逆矩阵$\mathbf{R}^{-1}$ $$ \mathbf{i} = \mathbf{R}^{-1} \mathbf{v}, $$ ### 多电子系统电子结构这里以 $H_3^+$ 为例，我们可以得到一个非常有用的系统。三个质子组成了一个等边三角形，并且共享了系统中的两个电子。假设我们可以用一个简单的哈密顿量 $\mathcal{H}$ 来描述这个系统，包含一项描述电子从一个位置跃迁到另一个位置，以及另一项用于描述双占据位置上两个电子之间的排斥： $$ \mathcal{H} = -t\sum_{i \neq j} a^\dagger_{i \sigma} a_{j \sigma} + U \sum n_{i\uparrow} n_{i\downarrow}, $$ 其中$a^\dagger_{i \sigma}$和$a_{j \sigma}$代表自旋向上$(\sigma = \uparrow)$ 和向下$(\sigma = \downarrow)$的产生和湮灭算符。 $n_{i\sigma} = a^\dagger_{i \sigma}a_{j \sigma}$是第 $i$ 层上的占据数。当用于描述高度关联电子体系时，这个哈密顿量称为Hubbard模型。其中的$t$和$U$可以从量子化学或者实验中得到。系统的Schro\"dinger方程为 $$ \mathcal{H}\ket{\Psi_k} = \mathcal{E}_k \ket{\Psi_k}. $$ 因为每个位置只有一个轨道，所以在单粒子表示下，两个电子总共有15个可能的状态。 $$ \begin{aligned}[l] \ket{\phi_1} = a^\dagger_{1\uparrow} a^\dagger_{1\downarrow} \ket{0}, \ket{\phi_2} = a^\dagger_{2\uparrow} a^\dagger_{2\downarrow} \ket{0}, \ket{\phi_3} = a^\dagger_{3\uparrow} a^\dagger_{3\downarrow} \ket{0}, \\ \ket{\phi_4} = a^\dagger_{1\uparrow} a^\dagger_{2\downarrow} \ket{0}, \ket{\phi_5} = a^\dagger_{1\uparrow} a^\dagger_{3\downarrow} \ket{0}, \ket{\phi_6} = a^\dagger_{2\uparrow} a^\dagger_{3\downarrow} \ket{0}, \\ \ket{\phi_7} = a^\dagger_{1\downarrow} a^\dagger_{2\uparrow} \ket{0}, \ket{\phi_8} = a^\dagger_{1\downarrow} a^\dagger_{3\uparrow} \ket{0}, \ket{\phi_9} = a^\dagger_{2\downarrow} a^\dagger_{3\uparrow} \ket{0}, \\ \ket{\phi_10} = a^\dagger_{1\uparrow} a^\dagger_{2\uparrow} \ket{0}, \ket{\phi_11} = a^\dagger_{1\uparrow} a^\dagger_{3\uparrow} \ket{0}, \ket{\phi_12} = a^\dagger_{2\uparrow} a^\dagger_{3\uparrow} \ket{0}, \\ \ket{\phi_13} = a^\dagger_{1\downarrow} a^\dagger_{2\downarrow} \ket{0}, \ket{\phi_14} = a^\dagger_{1\downarrow} a^\dagger_{3\downarrow} \ket{0}, \ket{\phi_15} = a^\dagger_{2\downarrow} a^\dagger_{3\downarrow} \ket{0}, \end{aligned} $$ 其中$\ket{0}$为真空态。此时，Schrodinger方程可以写成矩阵形式： $$ \mathbf{H \Psi}_k = \mathcal{E}_k \mathbf{\Psi}_k $$ 其中$H_{ij} = \bra{\phi_i}\mathcal{H}\ket{\phi_j}$，其中$i, j = 1, 2, \dots, 15$，$\Psi_{ki} = \braket{\phi_i|\Psi_k}$，最后也可以归结为一个久期方程： $$ |\mathbf{H} - \mathcal{E}I| = 0 $$ 其中$\mathbf{I}$ 为单位矩阵 $I_{ij} = \delta_{ij}$。这里的量子问题与分子的经典振动具有相似的数学形式。我们可以利用系统对称性对问题进行简化。总自旋以及总自旋的$z$分量与哈密顿量对易，因此，它们是好量子数。此时，我们可以将哈密顿矩阵越化为最小块大小为$2 \times 2$的块对角矩阵。当得到系统所有的本征值和本征矢量，就可以分析$H_3^+$的电子、光学、磁学性质。 ---- ## Basic matrix operations + $$ \mathbf{Ax} = \mathbf{b} \leftrightarrow \sum_{j=1}^nA_{ij}x_j = b_i, $$ + $$ \mathbf{C} = \mathbf{AB} \leftrightarrow C_{ij} = \sum_{k=1}^nA_{ik}B_{kj}, $$ + $$ \mathbf{AA}^{-1} = \mathbf{A}^{-1}\mathbf{A} = \mathbf{I} $$ + Derminant: $$ |\mathbf{A}| = \sum_{i=1}^n (-1)^{i+j}A_{ij}|\mathbf{R}_{ij}| $$ $|\mathbf{R}_{ij}|$ is the determinant of the residual matrix $\mathbf{R}_{ij}$ of $\mathbf{A}$ with its $i$th row and $j$th column removed. $C_{ij} = (-1)^{i+j}|\mathbf{R}_{ij}|$ is called a cofactor of $\mathbf{A}$. + Inverse: $$ A^{-1}_{ij} = \frac{C_{ij}}{|\mathbf{A}|}, $$ + If $|\mathbf{A}| \neq 0$, it is called a nonsingular matrix. + If $|\mathbf{A}| = 0$, it is called a singular matrix. + Trace of $\mathbf{A}$: sum of all its diagonal elements, $tr(\mathbf{A}) = \sum_{i=1}^n A_{ii}$. + Transpose: $$ A^T_{ij} = A_{ji} $$ + Orthogonality: $$ \mathbf{A}^T\mathbf{A} = \mathbf{AA}^T = \mathbf{I} $$ + Hermitian: complex conjugate $\mathbf{A}^\dagger \equiv \mathbf{A}^{T*}$ $$ \mathbf{A}^\dagger = \mathbf{A}^{-1} $$ + Add one row/column multiplied by a factor $\lambda$ to another row/column $$ A'_{ij} = A_{ij} + \lambda A_{kj}, $$ This operation preserves the determinant, $|\mathbf{A}'| = |\mathbf{A}|$ and can be represented by a matrix multiplication $$ \mathbf{A}' = \mathbf{M}\mathbf{A} $$ 其中，$\mathbf{M}$具有单位对角元，再加上一个非0非对角元 $M_{ik} = \lambda$。当调换任意两个行或者列，行列式值会变号。 $H^+_3$电子结构可以定义为 $$ \mathbf{Ax} = \lambda \mathbf{x} $$ 其中 $\mathbf{x}$ 和 $\lambda$ 为矩阵的本征矢量以及相应的本征值。如果将 $\mathbf{M}^{-1}\mathbf{A} = \mathbf{B}$ 作为问题的矩阵，$\mathbf{Ax} = \lambda \mathbf{Mx}$，将有 $\mathbf{Bx} = \lambda \mathbf{x}$。矩阵本征值问题可以看作是一个迭代求解的线性方程组问题。例如如果我们想求解上述方程的本征值和本征矢量，可以使用递归法 $$ \mathbf{Ax}^{(k+1)} = \lambda^{(k)} + \mathbf{x}^{(k)} $$ 其中$\lambda^{(k)}$和$\mathbf{x}^{(k)}$ 为第$k$次迭代中合适的本征值以及本征矢量，即，迭代法。当我们选择$\mathbf{x} = \mathbf{Sy}$时，由于 $$ \mathbf{By} = \lambda\mathbf{y} $$ 与 $$ \mathbf{By} = \lambda\mathbf{y} $$ 等价，矩阵本征值在相似变换 $$ \mathbf{B} = \mathbf{S}^{-1}\mathbf{A}\mathbf{S} $$ 下与矩阵$\mathbf{A}$本征值之一相同，其中$\mathbf{S}$矩阵非奇异，同时也意味着 $$ |\mathbf{B}| = |\mathbf{A}| = \prod_{i=1}^n \lambda_i $$ ## 线性方程系统对角线以上、以下矩阵元为0的矩阵成为上三角、下三角矩阵。求解矩阵最简单的方法就是```高斯消元法```，就是首先使用基本矩阵操作将系数矩阵变换成上、下三角矩阵，再通过向后、前替换求解方程的解。矩阵的逆以及行列式也可以通过这种方式求的。以$$ \mathbf{Ax} = \mathbf{b}$$ 为例。如果$|\mathbf{A}| \neq 0$以及$\mathbf{b} \neq 0$，系统就具有唯一解。将原来的矩阵标记为$\mathbf{A}^(0)$，$\mathbf{A}^(j)$为$j$次操作后的矩阵。$\mathbf{b}$也采取同样的标记。 ### 高斯消元高斯消元的基本思想是将原来的矩阵变换为具有相同解的上、下三角矩阵，这里以变换成上三角为例。同样也可以通过交换行与列的角色来得到下三角矩阵。在每次操作中，我们消掉不属于上三角的列中元素。最终的线性方程组为 $$ \mathbf{A}^{(n-1)}\mathbf{x} = \mathbf{b}^{(n-1)} $$ 其中$A_{ij}^{(n-1)} = 0~~\text{for}~~i>j$，这个过程很简单： :::{.callout-note} ### Gauss elimination - 首先给第一个方程乘以$-A_{i1}^{(0)}/A_{11}^{(0)}$，也就是矩阵的第一行以及$b_1^{(0)}$，并将其加到第$i$个方程对于所有的$i>1$； - 除第一行外，矩阵中每行的第一个元素都被消去，矩阵此时变成$\mathbf{A}^{(1)}$； - 第二个方程乘以$-A_{i2}^{(1)}/A_{22}^{(1)}$，并将其加到其他方程$i>2$； - 除第一、二行的第二个元素，其他行的第二个元素都别消掉了； - 对第三、四、$n-1$个方程重复该过程，该系数矩阵就变成了上三角矩阵 $\mathbf{A}^{(n-1)}$； ::: 上、下三角系数矩阵可以通过向后代入轻易求解。 <red>由于所有的对角线元都会被用到分母中，当任意一个分母为0时，该方法将不再适用。这可以通过行调换解决：将 $k$ 行与 $k$ 行以下 $k$ 列中最大元素 $A_{rk}$ 所在行$r$行进行调换</red>。该过程只是调换了各个方程的顺序，不改变方程的解。这里我们考虑部分调换情况，即只对特定列的剩余行进行调换。完全调换方法中会在所有剩余元素中搜索调换元。这里以部分调换为例。我们先从$|A^{(0)}_{1i}|$中搜索最大的元素，假设为$A^{(0)}_{k_1 1}$，之后将第一行与$k_1$行调换，从而消除第一行以外每行的第一个元素。之后对每行都进行类似操作就可以将$A$矩阵变换成上三角矩阵。物理上，在寻找枢轴元时我们不需要交换矩阵中的行。可以用一个指标来记录所有行中枢轴元的顺序。此外，对于$i>j$时，我们需要$A^{(j-1)}_{ij} / A^{(j-1)}_{k_j j}$比率来调整$\mathbf{b}^{(j-1)}$到$\mathbf{b}^j$，并存储到与$\mathbf{A}^{(n-1)}$相同二维矩阵中$A_{ij}^{(n-1)} = 0$的位置，即，对角线以下。当我们确定枢轴元时，为了得到公平的比较时，我们需要将每行的元素用所在行最大的元素进行重新标度。这个重新标度同样减少了潜在的舍入误差。 ```{python} {{< include python/gauss_elimination.py >}} ``` ### 矩阵行列式得到上三角或者下三角后，就很容易得到行列式。高斯消元只改变行列式的符号，不改变值，而符号可以通过行列的轮换次数确定。对于上三角矩阵，其行列式由对角线元素乘积给出。 ```{python} res = 1.0 for i in range(len(BB)): res = res * BB[i, i] print(res * sig) print(np.linalg.det(AA)) print(np.linalg.det(BB) * sig) ``` ### 线性方程组解类似的，通过高斯消元法将矩阵进行对角化的方法也可以用来求解线性方程组，解可以通过向后代入求的： $$ x_i = \frac{1}{A^{(n-1)}_{k_i i}} \left( b^{(n-1)}_{k_i} - \sum_{j=i+1}^n A^{(n-1)}_{k_i j} x_j \right) $$ 对于 $i = n-1, n - 2, \dots, 1$，出发点为 $x_n = b_{k_n}^{(n-1)} / A_{k_n n}^{(n-1)}$。这里使用 $k_i$ 代表 $i$ 列枢轴元素的行指标。容易证实，上述递归满足线性方程组 $\mathbf{A}^{(n-1)}\mathbf{x} = \mathbf{b}^{(n-1)}$。这里我们以惠斯通电桥中求解电流为例，其中 $r_s = r_1 = r_2 = r_x = r_a = 100~\Omega$ 以及 $v_0 = 200~\text{V}$。 ```{python} AA = np.array([[100.0, 100.0, 100.0], [-100.0, 300.0, -100.0], [-100.0, -100.0, 300.0]]) b = np.array([200.0, 0.0, 0.0]).reshape(3, 1) problem = np.concatenate((AA, b), axis = 1) print(problem.shape) print(problem) res, sig = gauss_elim(problem) A = res[:, 0:3] B = res[:, -1] x = np.zeros(len(AA)) x[-1] = B[-1] / A[-1, -1] for n in range(len(AA)-2, -1, -1): sm = 0.0 for i in range(n+1, len(AA)): sm += A[n, i] * x[i] x[n] = 1.0 / A[n, n] * (B[n] - sm) print(x) ``` ### 矩阵的逆 $\mathbf{A}$ 的逆矩阵 $\mathbf{A}^{-1}$ 可以通过类似方式求的： $$ A^{-1}_{ij} = x_{ij}, $$ 其中 $x_{ij}$ 为方程 $\mathbf{A}\mathbf{x}_j = \mathbf{b}$ 的解，且，$b_{ij} = \delta_{ij}$。如果将 $\mathbf{b}$ 展开成单位矩阵来求解线性方程，单位矩阵每列的解组成了 $\mathbf{A}^{-1}$ 的逆矩阵，此时只需要对上面程序进行稍加修改即可。 ```{python} def solveeq(AA, BB): problem = np.concatenate((AA, BB), axis = 1) res, sig = gauss_elim(problem) A = res[:, 0:3] B = res[:, -1] x = np.zeros(len(A)) x[-1] = B[-1] / A[-1, -1] for n in range(len(A)-2, -1, -1): sm = 0.0 for i in range(n+1, len(A)): sm += A[n, i] * x[i] x[n] = 1.0 / A[n, n] * (B[n] - sm) return x AA = np.array([[100.0, 100.0, 100.0], [-100.0, 300.0, -100.0], [-100.0, -100.0, 300.0]]) def inverse(AA): result = np.zeros(AA.shape) for i in range(len(AA)): BB = np.zeros(len(AA)).reshape(len(AA), 1) BB[i] = 1.0 result[:, i] = solveeq(AA, BB) return result result = inverse(AA) print(AA) print(result) print(np.dot(AA, result)) print(np.linalg.inv(AA)) ``` ### LU分解与高斯消元法相关的方法是LU分解，其将矩阵分解成一个上三角和一个下三角矩阵的乘积。如前所属，高斯消元法对应于$n-1$步矩阵操作并且可以表示为矩阵乘法 $$ \mathbf{A}^{(n-1)} = \mathbf{MA}, $$ {#eq-ama} 其中 $$ \mathbf{M} = \mathbf{M}^{(n-1)}\mathbf{M}^{(n-2)}\dots\mathbf{M}^{(1)}, $$ 其中每个 $\mathbf{M}^{(r)},~~\text{for}~~ r = 1, 2, \dots, n-1$，代表一步消元。容易验证 $\mathbf{M}$ 是下三角矩阵其对角元为 $-1$, $M_{k_ij} = -A_{k_i j}^{(j-1)} / A_{k_j j}^{(j-1)}~~~\text{for}~~~ i > j$，方程(@eq-ama)等价于 $$ \mathbf{A} = \mathbf{LU}. $$ 其中 $\mathbf{U} = \mathbf{A}^{(n-1)}$ 为上三角矩阵，$\mathbf{L} = \mathbf{M}^{-1}$ 为下三角矩阵。注意下三角矩阵的逆矩阵仍然为一个下三角矩阵。由于 $M_{ii} = -1$，$L_{ii} = 1$，$L_{ij} = -M_{ij} = A_{k_ij}^{(j-1)}/A_{k_jj}^{(j-1)} ~~~ \text{for} ~~~ i>j$，这就是之前用于部分轮换高斯消元时存储的比例。这个 $\mathbf{LU}$ 分解的方法将 $\mathbf{U}$ 的非 0 元素以及 $\mathbf{L}$ 的非0，非对角元素存储到了返回矩阵。 LU分解中选择 $L_{ii} = 1$ 的方法称为 Doolittle 分解。另外一个方法是选择 $U_{ii} = 1$；这称为Crout分解。Crout分解等价于用 $U_{ii}$ 将 $U_{ij}$ 缩放对于 $i<j$，以及给 $L_{ij}$ 乘以 $U_{jj}$ 对于 $i>j$（ $L_{ij}$ 来自于Doolitte分解）。一般情况下，$\mathbf{L}$ 和 $\mathbf{U}$ 的元素可以从二者乘积与 $\mathbf{A}$ 对比得到，结果可以写成一个迭代形式： $$ \begin{aligned}[l] L_{ij} = \frac{1}{U_{ij}}\left(A_{ij} - \sum_{k=1}^{j-1}L_{ik}U_{kj}\right), \\ U_{ij} = \frac{1}{L_{ij}}\left(A_{ij} - \sum_{k=1}^{i-1}L_{ik}U_{kj}\right), \end{aligned} $$ 其中 $L_{i1} = A_{i1}/U_{11}, U_{1j} = A_{1j} / L_{11}$ 为起始点。以此作为迭代求解，仍然需要考虑轮换的问题。