求助：用Crank-Nicolson法结合np.linalg.solve求解波动方程

一、线性方程组的直接求解 #

给定方程组（x、y为待求向量，c、d为已知向量，A、B为矩阵）：

1. x - y = c + d
2. y - A x = B d + c

消元求解步骤： #

• 从第一个方程解出 y = x - c - d，代入第二个方程：

  (x - c - d) - A x = B d + c

• 整理含x的项与常数项（I为单位矩阵）：

  (I - A)x = (B + I)d + 2c

• 若矩阵(I - A)可逆，直接得到x的解：

  x = (I - A)⁻¹[(B + I)d + 2c]

• 将x代回y = x - c - d，得到y的解：

  y = (I - A)⁻¹[(B + I)d + 2c] - c - d

二、Crank-Nicolson方法在波动方程中的应用关联 #

你熟悉该方法在一阶时间导数的含时薛定谔方程中的应用，波动方程是二阶时间导数，离散后会转化为上述形式的方程组，核心对应逻辑如下：

• 波动方程标准形式（以齐次为例）：∂²u/∂t² = ∇²u
• Crank-Nicolson方法对时间、空间做中心差分：
  • 时间二阶导数用三层差分：(uⁿ⁺¹ - 2uⁿ + uⁿ⁻¹)/Δt²
  • 空间拉普拉斯算子取tⁿ与tⁿ⁺¹时刻的平均：(∇²uⁿ⁺¹ + ∇²uⁿ)/2
• 整理后会得到关联uⁿ⁺¹（对应x）、uⁿ⁻¹（对应y）的线性方程组，其中c、d是包含uⁿ和源项的已知向量，A、B是空间离散后的组合矩阵（由拉普拉斯差分格式与时间步长参数构成）。

核心思路总结： #

• 波动方程的二阶时间导数离散后涉及三个时间层，Crank-Nicolson通过空间算子的时间平均，将其转化为双未知向量的线性方程组，本质和薛定谔方程的应用逻辑一致——隐式方法结合时间层平均，保证二阶精度且无条件稳定。
• 求解流程和薛定谔方程完全相通：将方程组消元为单未知向量的矩阵方程后，通过矩阵求逆或迭代法（如共轭梯度法）求解即可。

内容的提问来源于stack exchange，提问作者Guillermo Fuentes