9 Ordinary Differential Equations

  • 初值问题: 在某个时间点给定系统的状态, 求解系统在未来时间的演化
  • 边值问题: 在空间或时间的边界上给定系统的状态, 求解系统在整个区域内的行为

9.1 Euler Method

考虑一阶常微分方程的初值问题

\[\frac{\mathrm{d} y}{\mathrm{d} t}=f(t,y),\qquad a\leq t\leq b,\qquad y(a)=\alpha .\]

  • 在等间距网格

    \[t_i=a+i h,\qquad i=0,1,\cdots,N,\qquad h=\frac{b-a}{N}\]

    上, 数值解记为 \(w_i\approx y(t_i)\). 在 \((t_i,t_{i+1})\) 上, 对精确解作 Taylor 展开:

    \[\begin{aligned} y(t_{i+1}) &= y(t_i)+h y'(t_i)+\frac{h^2}{2}y''(\xi_i)\\ &= y(t_i)+h f(t_i,y(t_i))+\frac{h^2}{2}y''(\xi_i) \qquad \xi_i\in(t_i,t_{i+1}) \end{aligned}\]

  • Euler 法忽略余项, 并把精确值 \(y(t_i)\) 替换为数值近似 \(w_i\), 得到直接递归迭代

    \[w_0=\alpha,\qquad w_{i+1}=w_i+h f(t_i,w_i),\qquad i=0,1,\cdots,N-1.\]

9.2 Higher-Order Taylor Expansion

Euler 法是 Taylor 法的一阶截断. 若保留更高阶导数, 有

\[\begin{aligned} y(t_{i+1}) &=y(t_i)+h y'(t_i)+\frac{h^2}{2!}y''(t_i)+\cdots+\frac{h^n}{n!}y^{(n)}(t_i) +\frac{h^{n+1}}{(n+1)!}y^{(n+1)}(\xi_i)\\ &=y(t_i)+h T^{(n)}(t_i,y(t_i)) +\frac{h^{n+1}}{(n+1)!}f^{(n)}(\xi_i,y(\xi_i)) \end{aligned}\]

  • 数值格点中

    \[T^{(n)}(t_i,w_i) =f(t_i,w_i)+\frac{h}{2}f'(t_i,w_i)+\cdots+\frac{h^{n-1}}{n!}f^{(n-1)}(t_i,w_i)\]

    \(f',f'',\cdots\) 表示沿解轨道的全导数. 于是 \(n\) 阶 Taylor 迭代为

    \[w_0=\alpha,\qquad w_{i+1}=w_i+hT^{(n)}(t_i,w_i),\qquad i=0,1,\cdots,N-1\]

  • 局部截断误差

    \[\tau_{i+1}(h)=\frac{y(t_{i+1})-y(t_i)-hT^{(n)}(t_i,y(t_i))}{h}=\frac{h^n}{(n+1)!}f^{(n)}(\xi_i,y(\xi_i))\]

    按这里除以 \(h\) 的定义, 量级为 \(O(h^n)\); 未除以 \(h\) 的单步缺陷为 \(O(h^{n+1})\). 若导数充分光滑且 \(f\) 对 \(y\) 满足 Lipschitz 条件, 固定积分区间上的全局误差为 \(O(h^n)\).

Example 考虑初值问题

\[y'=y-t^2+1,\qquad 0\leq t\leq 2,\qquad y(0)=0.5.\]

精确解为 \(y(t)=(t+1)^2-\frac{1}{2}\mathrm{e}^{t}\).

Solution 取二阶 Taylor 法时(注意全导数)

\[T^{(2)}(t_i,w_i) =w_i-t_i^2+1+\frac{h}{2}(w_i-t_i^2+1-2t_i)\]

取步长 \(h=0.2\), \(t_i=0.2i\), 则

\[\begin{aligned} w_{i+1} &=w_i+0.2\left[\left(1+\frac{0.2}{2}\right)(w_i-0.04i^2+1)-0.04i\right]\\ &=1.22w_i-0.0088i^2-0.008i+0.22. \end{aligned}\]

高阶 Taylor 法可以显著减小截断误差, 但代价是需要计算 \(f\) 沿解轨道的高阶全导数. 因此实际上常不如 Runge-Kutta 类方法方便.

9.3 Runge-Kutta Methods

Runge-Kutta 方法的思想是用若干次函数值计算来代替 Taylor 法中的显式高阶导数.

9.3.1 Midpoint Method

从二阶 Taylor 展开出发

\[T^{(2)}(t,y)=f(t,y)+\frac{h}{2}f'(t,y),\]

其中沿解曲线的全导数为

\[f'(t,y)=\frac{\mathrm{d}f}{\mathrm{d}t}(t,y) =\frac{\partial f}{\partial t}(t,y)+\frac{\partial f}{\partial y}(t,y)y'(t) =\frac{\partial f}{\partial t}(t,y)+\frac{\partial f}{\partial y}(t,y)f(t,y).\]

  • 若用一个移位函数值 \(a_1 f(t+\alpha_1,y+\beta_1)\) 近似 \(T^{(2)}\), 由一阶展开

    \[a_1 f(t+\alpha_1,y+\beta_1) \approx a_1 f(t,y)+a_1\alpha_1\frac{\partial f}{\partial t}(t,y) +a_1\beta_1\frac{\partial f}{\partial y}(t,y)\]

    与 \(T^{(2)}\) 的系数比较可取

    \[a_1=1,\qquad \alpha_1=\frac{h}{2},\qquad \beta_1=\frac{h}{2}f(t,y).\]

    以对其到 1 阶全导数.

  • 忽略余项时得到二阶中点法

    \[w_{i+1}=w_i+h f\left(t_i+\frac{h}{2},w_i+\frac{h}{2}f(t_i,w_i)\right).\]

9.3.2 Improved Euler Method

另一个二阶形式是改进 Euler 法. 由

\[y''(t)=\frac{\mathrm{d}}{\mathrm{d}t}f(t,y) =f_t(t,y)+f_y(t,y)y' =f_t+f_y f\]

可得

\[y(t+h)=y(t)+h f+\frac{h^2}{2}(f_t+f_y f)+O(h^3).\]

  • 尝试取两个函数值的线性组合作为增率函数:

    \[\begin{aligned} \Phi &\equiv a_1 f(t,y)+a_2 f(t+\alpha h,y+\beta h f)\\ &=(a_1+a_2)f+a_2\alpha hf_t+a_2\beta hf_yf+\mathcal{O}(h^2) \end{aligned}\]

    比较系数, 取 \(a_1=a_2=1/2,\alpha=\beta=1\). 于是

    \[\Phi=\frac{1}{2}f(t,y)+\frac{1}{2}f(t+h,y+h f).\]

  • 结果上, 改进 Euler 法使用初始斜率和 Euler 预测端点斜率的平均:

    \[w_{i+1}=w_i+\frac{h}{2}\left[f(t_i,w_i)+f(t_{i+1},w_i+h f(t_i,w_i))\right].\]

9.3.3 Third-Order Runge-Kutta Method

更高阶 Runge-Kutta 法继续把 Taylor 目标中的导数信息换成多次函数值. 例如三阶目标为

\[T^{(3)}(t,y)=f(t,y)+\frac{h}{2}f'(t,y)+\frac{h^2}{6}f''(t,y).\]

经过推导可得到三阶公式

\[\begin{aligned} w_0&=\alpha,\\ k_1&=h f(t_i,w_i),\\ k_2&=h f\left(t_i+\frac{h}{3},w_i+\frac{k_1}{3}\right),\\ k_3&=h f\left(t_i+\frac{2h}{3},w_i+\frac{2k_2}{3}\right),\\ w_{i+1}&=w_i+\frac{1}{4}(k_1+3k_3), \end{aligned}\]

其中 \(i=0,1,\cdots,N-1\).

9.3.4 Fourth-Order Runge-Kutta Method

经典四阶 Runge-Kutta 方法, 也就是 RK4, 为

\[\begin{aligned} w_0&=\alpha,\\ k_1&=h f(t_i,w_i),\\ k_2&=h f\left(t_i+\frac{h}{2},w_i+\frac{k_1}{2}\right),\\ k_3&=h f\left(t_i+\frac{h}{2},w_i+\frac{k_2}{2}\right),\\ k_4&=h f(t_{i+1},w_i+k_3),\\ w_{i+1}&=w_i+\frac{1}{6}(k_1+2k_2+2k_3+k_4). \end{aligned}\]

RK4 在一般初值问题中使用广泛, 因为它不需要显式求高阶导数, 同时具有四阶精度. 对向量初值问题, 将 \(w_j,f,k_s\) 换成向量即可使用相同权重; 每个 stage 必须使用同一中间状态计算全部分量, 不能逐分量更新后混用新旧状态.

9.4 Predictor-Corrector Adams Methods *

多步法利用多个历史点来外推或插值 \(f(t,y)\), 再对一步区间上的积分进行近似. 与单步 RK(Runge-Kutta) 法相比, 它在每一步可以复用过去已经计算过的斜率, 但需要若干初始点启动.

9.4.1 Adams-Moulton Fourth-Order Implicit Method

三步四阶隐式 Adams-Moulton 格式为

\[\begin{aligned} w_0&=\alpha,\qquad w_1=\alpha_1,\qquad w_2=\alpha_2,\\ F_i&=f(t_i,w_i),\\ w_{i+1}&=w_i+\frac{h}{24}(9F_{i+1}+19F_i-5F_{i-1}+F_{i-2}). \end{aligned}\]

这里 \(w_{i+1}\) 同时出现在左边和右边, 因而是 implicit method.

9.4.2 Adams-Bashforth Third-Order Explicit Method

相应的三步显式 Adams-Bashforth 格式为

\[\begin{aligned} w_0&=\alpha,\qquad w_1=\alpha_1,\qquad w_2=\alpha_2,\\ F_i&=f(t_i,w_i),\\ w_{i+1} &=w_i+\frac{h}{12}(23F_i-16F_{i-1}+5F_{i-2}). \end{aligned}\]

9.4.3 Adams-Bashforth Fourth-Order Explicit Method

四步显式 Adams-Bashforth 格式为

\[\begin{aligned} w_0&=\alpha,\qquad w_1=\alpha_1,\qquad w_2=\alpha_2,\qquad w_3=\alpha_3,\\ F_i&=f(t_i,w_i),\\ w_{i+1} &=w_i+\frac{h}{24}(55F_i-59F_{i-1}+37F_{i-2}-9F_{i-3}). \end{aligned}\]

实际计算时, 可以先用四阶 Runge-Kutta 方法计算 \(w_0,w_1,w_2,w_3\). 之后用四步 Adams-Bashforth 作 predictor:

\[w_{4p}=w_3+\frac{h}{24}(55F_3-59F_2+37F_1-9F_0).\]

再用三步隐式 Adams-Moulton 作 corrector, 把未知的 \(w_4\) 先由 predictor \(w_{4p}\) 代入右端:

\[F_{4p}=f(t_4,w_{4p}),\qquad w_4=w_3+\frac{h}{24}(9F_{4p}+19F_3-5F_2+F_1).\]

下一步使用修正值重新计算 \(F_4=f(t_4,w_4)\), 再继续推进. 这里只作一次校正, 因而整个 predictor-corrector 递推是显式的, 并没有严格求解隐式 Adams-Moulton 方程. 在 \(f\) 光滑、启动值具有四阶精度的条件下, 它仍保持四阶精度, 但稳定域与完全求解的隐式格式不同.

9.5 Extended Fourth-Order Runge-Kutta Method *

常微分方程组的初值问题写成

\[\frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t} =\mathbf{f}(t,\mathbf{u}), \qquad \mathbf{u}(a)=\boldsymbol{\alpha}, \qquad \mathbf{u}=(u_1,\cdots,u_m)^{\mathsf{T}}, \qquad a\leq t\leq b.\]

单变量 RK4 的四个 stage 和权重可直接推广到向量. 设

\[\mathbf{w}_j=(w_{1,j},\cdots,w_{m,j})^{\mathsf{T}} \approx \mathbf{u}(t_j), \qquad \mathbf{w}_0=\boldsymbol{\alpha}.\]

按 stage 定义

\[\begin{aligned} \mathbf{k}_1&=h\mathbf{f}(t_j,\mathbf{w}_j),\\ \mathbf{k}_2&=h\mathbf{f}\left(t_j+\frac{h}{2},\mathbf{w}_j+\frac{\mathbf{k}_1}{2}\right),\\ \mathbf{k}_3&=h\mathbf{f}\left(t_j+\frac{h}{2},\mathbf{w}_j+\frac{\mathbf{k}_2}{2}\right),\\ \mathbf{k}_4&=h\mathbf{f}(t_j+h,\mathbf{w}_j+\mathbf{k}_3). \end{aligned}\]

递推公式为

\[\mathbf{w}_{j+1} =\mathbf{w}_j+\frac{1}{6}(\mathbf{k}_1+2\mathbf{k}_2+2\mathbf{k}_3+\mathbf{k}_4).\]

第 \(\ell\) 个分量就是原来的 \(k_{s,\ell}\) 公式. 由于 \(\mathbf{k}_2\), \(\mathbf{k}_3\), \(\mathbf{k}_4\) 都依赖前一 stage 中所有变量的增量, 方程组必须按 stage 整体推进.

9.6 Higher-Order Ordinary Differential Equations

高阶方程

\[y^{(m)}(t)=f(t,y,y',\cdots,y^{(m-1)}),\qquad a\leq t\leq b,\]

初值

\[y(a)=\alpha_1,\qquad y'(a)=\alpha_2,\qquad \cdots\qquad y^{(m-1)}(a)=\alpha_m\]

也可化为一阶.

  • 令 \[u_1(t)=y(t),\qquad u_2(t)=y'(t),\qquad \cdots,\qquad u_m(t)=y^{(m-1)}(t),\] 则 \[\begin{cases} u_1'=u_2,\\ u_2'=u_3,\\ \cdots\\ u_{m-1}'=u_m,\\ u_m'=f(t,u_1,u_2,\cdots,u_m). \end{cases}\] 相应初值为 \[u_1(a)=\alpha_1,\qquad u_2(a)=\alpha_2,\qquad \cdots,\qquad u_m(a)=\alpha_m.\]

Example 将二阶方程化为一阶系统

\[y''-2y'+2y=\mathrm{e}^{2t}\sin t,\qquad 0\leq t\leq 1,\qquad y(0)=-0.4,\qquad y'(0)=-0.6\]

Solution 取

\[u_1(t)=y(t),\qquad u_2(t)=y'(t),\]

得到

\[\begin{aligned} u_1'(t)&=u_2(t),\\ u_2'(t)&=\mathrm{e}^{2t}\sin t-2u_1(t)+2u_2(t), \end{aligned} \qquad u_1(0)=-0.4,\quad u_2(0)=-0.6.\]

试用 \(h=0.1\) 的四阶 Runge-Kutta 求第一步, 初始值为 \(w_{1,0}=-0.4\), \(w_{2,0}=-0.6\), \(t_0=0\). stage 增量为

\[\mathbf{F}(t,\mathbf{u}) = \begin{pmatrix} u_2\\ \mathrm{e}^{2t}\sin t-2u_1+2u_2 \end{pmatrix}.\]

代入 RK4 stage 后得到

\[\begin{aligned} \mathbf{k}_1&=(-0.06,\,-0.04)^{\mathsf{T}},\\ \mathbf{k}_2&=(-0.062,\,-0.03247644756)^{\mathsf{T}},\\ \mathbf{k}_3&=(-0.06162382238,\,-0.03152409232)^{\mathsf{T}},\\ \mathbf{k}_4&=(-0.06315240923,\,-0.02178637294)^{\mathsf{T}}. \end{aligned}\]

因此

\[\mathbf{w}_1 =\mathbf{w}_0+\frac{1}{6}(\mathbf{k}_1+2\mathbf{k}_2+2\mathbf{k}_3+\mathbf{k}_4) =(-0.4617333423,\,-0.6316312421)^{\mathsf{T}}.\]

9.7 Stiff Differential Equations

stiff equations 包含衰减尺度相差很大的模式, 因而真实解中快速衰减的分量可能很小, 但显式格式仍会被这些快模态的稳定性限制住. 此时显式方法需要极小步长, 否则很容易发散不收敛. 隐式方法的稳定性更适合处理此类问题.

对测试方程 \(y'=\lambda y\), 令 \(z=h\lambda\), 显式 Euler 的放大因子为 \(R(z)=1+z\), 稳定条件为 \(|1+z|\le1\); 隐式 Euler 为 \(R(z)=(1-z)^{-1}\), 在整个左半平面稳定且 \(z\to-\infty\) 时趋于零. 因而负实快模态会严格限制显式 Euler 的步长, 而隐式 Euler 可以将其强烈衰减.

9.7.1 Implicit Trapezoidal Method

隐式梯形法使用梯形公式近似:

\[w_0=\alpha,\qquad w_{j+1}=w_j+\frac{h}{2}\left[f(t_{j+1},w_{j+1})+f(t_j,w_j)\right].\]

因为 \(w_{j+1}\) 出现在右端, 每一步需要求解

\[F(w)=w-w_j-\frac{h}{2}\left[f(t_{j+1},w)+f(t_j,w_j)\right]=0.\]

Newton method 的迭代写为

\[w_{j+1}^{(k)} =w_{j+1}^{(k-1)}-\frac{F(w_{j+1}^{(k-1)})}{F'(w_{j+1}^{(k-1)})}, \qquad F'(w)=1-\frac{h}{2}f_y(t_{j+1},w).\]

对测试方程 \(y'=\lambda y\), 梯形法的放大因子为 \(R(z)=(1+z/2)/(1-z/2)\), 其中 \(z=h\lambda\). 当 \(\operatorname{Re}z\le0\) 时有 \(|R(z)|\le1\), 即 A 稳定; 但 \(z\to-\infty\) 时 \(R(z)\to-1\), 对极快衰减模式可能产生数值振荡, 不具备 L 稳定性.

9.7.2 Gear Method *

Gear 方法使用后向差分公式 BDF(backward differentiation formula). 以下取固定步长 \(h\), 并把当前值的系数归一为 1. 一般写成

\[\mathbf{y}_{n+k}-h\beta_k\mathbf{f}(x_{n+k},\mathbf{y}_{n+k})=\mathbf{g}_{n+k}, \qquad \mathbf{g}_{n+k}=-\sum_{j=0}^{k-1}\alpha_j\mathbf{y}_{n+j}.\]

固定步长 BDF 在 1–6 阶满足零稳定性; 其中一阶和二阶为 A 稳定, 三至六阶为 \(A(\alpha)\) 稳定, 稳定域包含负实轴周围的扇形区域, 但不包含整个左半平面, 不能对所有左半平面特征值任意取步长. 更高阶 BDF 不满足零稳定性.1

二阶和三阶例子为

\[\begin{aligned} y_{n+2}&=\frac{4y_{n+1}-y_n+2hf_{n+2}}{3},\\ y_{n+3}&=\frac{18y_{n+2}-9y_{n+1}+2y_n+6hf_{n+3}}{11}. \end{aligned}\]

对向量问题, 当前步未知量可用 Newton method 求解:

\[\mathbf{y}_{n+k}^{[s+1]} =\mathbf{y}_{n+k}^{[s]} -\left[\mathbf{I}-h\beta_k\frac{\partial \mathbf{f}}{\partial \mathbf{y}}(x_{n+k},\mathbf{y}_{n+k}^{[s]})\right]^{-1} \left[\mathbf{y}_{n+k}^{[s]}-h\beta_k\mathbf{f}(x_{n+k},\mathbf{y}_{n+k}^{[s]})-\mathbf{g}_{n+k}\right].\]

9.8 Stability of Explicit and Implicit Euler Methods

Example 测试方程

\[y'=-cy,\qquad c>0 \implies y(t)=y_0\mathrm{e}^{-ct}.\]

Solution 显式差分给出

\[y_{n+1}=y_n+h y_n'=(1-ch)y_n \implies y_n=(1-ch)^n y_0.\]

若要求固定步长下数值解随 \(n\to\infty\) 衰减, 必须有 \(0\lt ch\lt 2\). 这是绝对稳定条件; 对固定时间区间逼近精确解的收敛则指 \(h\to0\). \(c\) 很大时, 显式 Euler 法要求步长 \(h\) 很小, 否则会不稳定.

而隐式差分给出

\[y_{n+1}=y_n+h y_{n+1}'=y_n-chy_{n+1} \implies y_{n+1}=\frac{y_n}{1+ch},\qquad y_n=\frac{y_0}{(1+ch)^n}.\]

因此隐式 Euler 法对该测试方程无条件稳定, 更适合刚性方程.

Example 上例推广到线性系统 \(\mathbf{y}'=-\mathbf{C}\mathbf{y}\).

Solution 先求严格解: 若 \(\mathbf{T}^{-1}\mathbf{C}\mathbf{T}=\operatorname{diag}(\lambda_0,\cdots,\lambda_{N-1})\), 令

\[\mathbf{z}=\mathbf{T}^{-1}\mathbf{y} \implies \mathbf{z}'=-\operatorname{diag}(\lambda_0,\cdots,\lambda_{N-1})\mathbf{z}\]

解为

\[\mathbf{y}(x)=\mathbf{T}\operatorname{diag}\left(\mathrm{e}^{-\lambda_0 x},\cdots,\mathrm{e}^{-\lambda_{N-1}x}\right)\mathbf{T}^{-1}\mathbf{y}_0.\]

当所有 \(\lambda_i\) 的实部为正时, 连续系统稳定. 刚性系统的特征值实部相差很大, 因而显式方法的步长会被最大特征值限制.

  • 显式差分满足

    \[\mathbf{y}_{n+1}=(\mathbf{I}-h\mathbf{C})\mathbf{y}_n,\qquad \mathbf{y}_n=(\mathbf{I}-h\mathbf{C})^n\mathbf{y}_0.\]

    若沿特征向量展开

    \[\mathbf{y}_0=\sum_{i=0}^{N-1}\alpha_i\boldsymbol{\xi}_i \implies \mathbf{y}_n=\sum_{i=0}^{N-1}\alpha_i(1-h\lambda_i)^n\boldsymbol{\xi}_i.\]

    稳定条件为 \(|1-h\lambda_i|\lt 1\); 若特征值均为正实数, 则化为

    \[|1-h\lambda_i|\lt 1 \implies h\lt \frac{2}{\lambda_{\max}}.\]

  • 隐式差分则给出

    \[\mathbf{y}_{n+1}=(\mathbf{I}+h\mathbf{C})^{-1}\mathbf{y}_n, \qquad |1+h\lambda_i|^{-1}\lt 1,\qquad i=0,\cdots,N-1.\]

    没有最大特征值步长限制.

9.9 Linear Shooting Method

考虑二阶线性边值问题

\[y''=p(x)y'+q(x)y+r(x),\qquad a\leq x\leq b, \qquad y(a)=\alpha,\qquad y(b)=\beta.\]

打靶法把未知初始斜率 \(y'(a)\) 视为待定参数 \(t\), 通过调整 \(t\) 使解命中右端边界, 从而把边值问题转化为初值问题.

在线性情形, 可用叠加原理将原解拆成两部分:

  1. 非齐次特解 \(y_1\): 满足原方程, 取零初始斜率

    \[y_1''=p(x)y_1'+q(x)y_1+r(x),\qquad y_1(a)=\alpha,\quad y_1'(a)=0.\]

  2. 齐次补解 \(y_2\): 满足齐次方程, 取单位初始斜率

    \[y_2''=p(x)y_2'+q(x)y_2,\qquad y_2(a)=0,\quad y_2'(a)=1.\]

则任意组合

\[y(x)=y_1(x)+t\cdot y_2(x)\]

都满足原方程和左边界 \(y(a)=\alpha\) (因为 \(y_2(a)=0\)), 且初始斜率恰好是 \(y'(a)=t\). 代入右边界:

\[y(b)=y_1(b)+t\,y_2(b)=\beta \;\Longrightarrow\; t=\frac{\beta-y_1(b)}{y_2(b)}.\]

若 \(y_2(b)\ne0\), 原问题的唯一解为

\[y(x)=y_1(x)+\frac{\beta-y_1(b)}{y_2(b)}\,y_2(x).\]

若 \(y_2(b)=0\), 则 \(\beta\ne y_1(b)\) 时无解, \(\beta=y_1(b)\) 时初始斜率不能由右边界唯一确定, 存在一族解. 这些结论假设两个初值问题在整个区间上存在唯一解.

实际计算时分别算出 \(y_1(x)\) 和 \(y_2(x)\), 再按上式混合. 这就是线性打靶法.

9.10 Nonlinear Shooting Method *

考虑二阶非线性边值问题

\[y''=f(x,y,y'),\qquad a\leq x\leq b, \qquad y(a)=\alpha,\qquad y(b)=\beta.\]

由于 \(x=a\) 处的导数未知, 设初始斜率 \(y'(a)=t\) 为待定参数, 转化为初值问题

\[y''=f(x,y,y'),\qquad y(a)=\alpha,\qquad y'(a)=t.\]

目标是调节 \(t\), 使数值解出的 \(y(b,t)\) 满足

\[F(t)=y(b,t)-\beta=0.\]

  • 割线法:

    \[t_k=t_{k-1} -\frac{F(t_{k-1})(t_{k-1}-t_{k-2})} {F(t_{k-1})-F(t_{k-2})}, \qquad k=2,3,\cdots.\]

  • Newton method

    \[t_k=t_{k-1}-\frac{F(t_{k-1})}{F'(t_{k-1})}\]

    但需要知道解对初始斜率的灵敏度. 对 \(y''(x,t)=f(x,y(x,t),y'(x,t))\) 关于参数 \(t\) 求导, 得

    \[\frac{\partial y''}{\partial t} =\frac{\partial f}{\partial y}\frac{\partial y}{\partial t} +\frac{\partial f}{\partial y'}\frac{\partial y'}{\partial t}.\]

    令

    \[z(x,t)=\frac{\partial y}{\partial t}(x,t), \qquad F'(t)=z(b,t),\]

    则灵敏度方程为

    \[z''(x,t)=f_y(x,y,y')z(x,t)+f_{y'}(x,y,y')z'(x,t), \qquad a\leq x\leq b.\]

    由于 \(y(a,t)=\alpha\) 不依赖 \(t\), 而 \(y'(a,t)=t\), 初值为

    \[z(a,t)=0,\qquad z'(a,t)=1.\]

    于是 Newton 递推可写为

    \[t_k=t_{k-1}-\frac{F(t_{k-1})}{z(b,t_{k-1})}.\]

    在每个试探 \(t\) 下, 同时求解 \(y(b,t)\) 和 \(z(b,t)\), 再用它们修正 \(t\). 因此 shooting method 把边值问题转化为一系列初值问题.

9.11 Numerov Method *

Numerov 方法针对

\[y''(x)=f(x)y(x)\]

这种不含一阶导数的二阶方程.

  • 从 Taylor 展开到六阶开始:

    \[y(x+h)+y(x-h)=2y(x)+h^2y''(x)+\frac{h^4}{12}y^{(4)}(x)+\frac{h^6}{360}y^{(6)}(x)+O(h^8).\]

  • 引入辅助函数

    \[w(x)=\left[1-\frac{h^2}{12}f(x)\right]y(x).\]

    用中心差分近似

    \[y^{(4)}(x)\approx\frac{y''(x+h)-2y''(x)+y''(x-h)}{h^2},\]

    可得 Numerov 递推

    \[w(x+h)+w(x-h)-2w(x)=h^2f(x)y(x)+O(h^6).\]

  • 起点公式为

    \[2h y'(0)=\left[1-\frac{h^2}{6}f(h)\right]y(h)-\left[1-\frac{h^2}{6}f(-h)\right]y(-h)+O(h^5).\]

    递推在起点给出

    \[\left[1-\frac{h^2}{12}f(h)\right]y(h)+\left[1-\frac{h^2}{12}f(-h)\right]y(-h) =\left[2+\frac{5h^2}{6}f(0)\right]y(0)+O(h^6).\]

    将该式乘以 \(1-h^2f(-h)/6\), 将起点导数式乘以 \(1-h^2f(-h)/12\), 再相加以消去 \(y(-h)\). 假设 \(f\) 在起点两侧充分光滑且下式的 \(B\ne0\), 得 \(y(h)=A/B+O(h^5)\), 其中

    \[\begin{aligned} A&=\left[2+\frac{5h^2}{6}f(0)\right]\left[1-\frac{h^2}{6}f(-h)\right]y(0)+2h y'(0)\left[1-\frac{h^2}{12}f(-h)\right],\\ B&=\left[1-\frac{h^2}{12}f(h)\right]\left[1-\frac{h^2}{6}f(-h)\right]+\left[1-\frac{h^2}{12}f(-h)\right]\left[1-\frac{h^2}{6}f(h)\right]. \end{aligned}\]

Numerov 方法也可用于

\[y''(x)+k(x)y(x)=F(x),\qquad a\leq x\leq b, \qquad y(a)=\alpha,\qquad y(b)=\beta.\]
  • 二阶导数可写作 \[y''(x)=\frac{y(x+h)+y(x-h)-2y(x)}{h^2}-\frac{h^2}{12}y^{(4)}(x)+O(h^4).\] 由 \(y''(x)=-k(x)y(x)+F(x)\), 得 \[y^{(4)}(x)=-\frac{\mathrm{d}^2}{\mathrm{d}x^2}\left(ky-F\right).\] 中心差分近似给出 \[g_i=k_i y_i-F_i,\qquad \frac{\mathrm{d}^2}{\mathrm{d}x^2}(ky-F)_i \approx\frac{g_{i+1}-2g_i+g_{i-1}}{h^2}.\] 代入并整理, 得递推公式 \[\left(1+\frac{h^2}{12}k_{i+1}\right)y_{i+1} -2\left(1-\frac{5h^2}{12}k_i\right)y_i +\left(1+\frac{h^2}{12}k_{i-1}\right)y_{i-1} =\frac{h^2}{12}(F_{i+1}+10F_i+F_{i-1})+O(h^6).\]

9.12 Green Function Method *

Green 函数法把线性边值问题 \(L_xu(x)=f(x)\) 写成线性响应形式.

若 \(G(x,s)\) 满足 \(L_xG(x,s)=\delta(x-s)\) 并满足同类齐次边界条件, 则解可以写为

\[u(x)=\int G(x,s)f(s)\,\mathrm{d}s.\]

Example

\[-u''(x)=f(x),\qquad x\in(0,1),\qquad u(0)=u(1)=0\]

Solution 作不定积分得到

\[\begin{aligned} u(x) &= c_1+c_2x-\int_0^x\left(\int_0^y f(z)\,\mathrm{d}z\right)\mathrm{d}y\\ &= c_1+c_2x-\int_0^x (x-z)f(z)\,\mathrm{d}z \end{aligned}\]

由边条件得

\[c_1=0, \qquad c_2=\int_0^1 (1-y)f(y)\,\mathrm{d}y.\]

于是

\[\begin{aligned} u(x) &= x\int_0^1(1-y)f(y)\,\mathrm{d}y-\int_0^x(x-y)f(y)\,\mathrm{d}y\\ &= \int_0^1G(x,y)f(y)\,\mathrm{d}y, \qquad G(x,y)= \begin{cases} y(1-x), & 0\leq y\leq x\\ x(1-y), & x\leq y\leq 1 \end{cases} \end{aligned}\]

  • \(G\) 受边条件

    \[G(0,y)=G(1,y)=G(x,0)=G(x,1)=0\]

    约束, 且在 \(x=y\) 处连续. 在 \(x=y\) 处导数有跳跃

    \[\left.\frac{\mathrm{d}G}{\mathrm{d}x}\right|_{x=y+\varepsilon} -\left.\frac{\mathrm{d}G}{\mathrm{d}x}\right|_{x=y-\varepsilon} =-1\]

  • 它可以用齐次方程 \(-u''=0\) 的两个解重新表达:

    1. \(u_\lt (x)=x\) 从左端出发, 满足左边界条件 \(u_\lt (0)=0\).
    2. \(u_>(x)=1-x\) 从右端出发, 满足右边界条件 \(u_>(1)=0\).

    当 \(y\leq x\) 时, \(G=y(1-x)=u_\lt (y)\,u_>(x)\); 当 \(x\leq y\) 时, \(G=x(1-y)=u_\lt (x)\,u_>(y)\). 归一化因子来自 \(L_xG(x,y)=\delta(x-y)\). 由 \(a_2=-1\)

    \[G(x,y)=\frac{1}{a_2(y)W(y)}\begin{cases} u_\lt (y)\,u_>(x), & y\leq x\\ u_\lt (x)\,u_>(y), & x\leq y \end{cases}\]

    其中 Wronskian 为

    \[W(y)=u_\lt (y)u_>'(y)-u_\lt '(y)u_>(y) =y\cdot(-1)-1\cdot(1-y)=-1.\]

    Green 函数由两个分别满足左右边界条件的齐次解拼成, 在源点 \(x=y\) 处通过 Wronskian 匹配. 并受导数跳跃条件约束.

在数值计算中, 若齐次方程不含一阶导数项, 可用 Numerov 分别从左向右积分得到 \(u_\lt \), 从右向左积分得到 \(u_>\). 若含一阶导数项, 须先通过变量变换消去该项, 或直接用适用的一阶系统积分法. 在齐次边值问题没有非零解的条件下, 计算 Wronskian 后在网格上构造 \(G(x_i,x_j)\), 最后用求积近似

\[u(x_i)\approx \sum_j w_jG(x_i,x_j)f(x_j).\]

一般的二阶线性边值问题 考虑

\[\frac{\mathrm{d}^2y}{\mathrm{d}x^2} +a_1(x)\frac{\mathrm{d}y}{\mathrm{d}x} +a_0(x)y=f(x), \qquad B_a[y]=0,\qquad B_b[y]=0.\]

其中 \(B_a,B_b\) 表示左右端的线性边界条件

  • Dirichlet 边界条件对应 \(B_a[y]=y(a)\), \(B_b[y]=y(b)\);
  • Robin 边界条件可写成 \(B_a[y]=\alpha_a y(a)+\beta_a y'(a)\), \(B_b[y]=\alpha_b y(b)+\beta_b y'(b)\).

相应齐次方程

\[\frac{\mathrm{d}^2y}{\mathrm{d}x^2} +a_1(x)\frac{\mathrm{d}y}{\mathrm{d}x} +a_0(x)y=0\]

设对应的齐次边值问题没有非零解. 取两个齐次解 \(y_1,y_2\), 使 \(y_1\) 满足 \(a\) 端边界条件, \(y_2\) 满足 \(b\) 端边界条件. 由常数变易法(variation of parameters)可写为

\[y_p(x)=c_1(x)y_1(x)+c_2(x)y_2(x),\]

可得

\[y(x)= \int_a^x \frac{y_1(s)y_2(x)f(s)}{W(y_1,y_2)(s)}\,\mathrm{d}s + \int_x^b \frac{y_1(x)y_2(s)f(s)}{W(y_1,y_2)(s)}\,\mathrm{d}s.\]

因此进一步给出解的叠加表示为

\[y(x)=\int_a^b G(x,s)f(s)\,\mathrm{d}s, \qquad G(x,s)= \begin{cases} \displaystyle \frac{y_1(s)y_2(x)}{W(y_1,y_2)(s)}, & a\leq s\leq x\leq b,\\[8pt] \displaystyle \frac{y_1(x)y_2(s)}{W(y_1,y_2)(s)}, & a\leq x\leq s\leq b, \end{cases}\]

Green 函数 \(G\) 是线性响应函数, 与外场 \(f\) 无关; 源项只通过最后的积分进入解.

9.13 Finite-Difference Method for ODEs

考虑边值问题

\[y''(x)+p(x)y'(x)+q(x)y(x)=r(x),\qquad a\leq x\leq b, \qquad y(a)=\alpha,\qquad y(b)=\beta.\]

  • 取 \(h=\frac{b-a}{N},x_n=a+nh\), 记

    \[y_n=y(x_n),\qquad p_n=p(x_n),\qquad q_n=q(x_n),\qquad r_n=r(x_n).\]

  • 中心差分公式为

    \[y'(x_n)=\frac{y_{n+1}-y_{n-1}}{2h}+O(h^2), \qquad y''(x_n)=\frac{y_{n+1}-2y_n+y_{n-1}}{h^2}+O(h^2).\]

    代入微分方程并整理为

    \[\ell_i y_{i-1}+d_i y_i+u_i y_{i+1}=r_i,\qquad i=1,\cdots,N-1.\]

    其中

    \[\ell_i=\frac{1}{h^2}-\frac{p_i}{2h},\qquad d_i=-\frac{2}{h^2}+q_i,\qquad u_i=\frac{1}{h^2}+\frac{p_i}{2h}.\]

  • 边界条件给出 \(y_0=\alpha,y_N=\beta\). 剩余未知量 \(\mathbf{y}=(y_1,\cdots,y_{N-1})^{\mathsf{T}}\) 化为三对角线性方程组

    \[\begin{pmatrix} d_1&u_1&0&\cdots&0\\ \ell_2&d_2&u_2&\cdots&0\\ 0&\ddots&\ddots&\ddots&\vdots\\ \vdots&\ddots&\ell_{N-2}&d_{N-2}&u_{N-2}\\ 0&\cdots&0&\ell_{N-1}&d_{N-1} \end{pmatrix} \begin{pmatrix} y_1\\y_2\\\vdots\\y_{N-2}\\y_{N-1} \end{pmatrix} = \begin{pmatrix} r_1-\ell_1\alpha\\r_2\\\vdots\\r_{N-2}\\r_{N-1}-u_{N-1}\beta \end{pmatrix}\]

    的求解问题.

10 Partial Differential Equations

10.1 Finite-Difference PDEs

偏微分方程包含多个独立变量. Poisson 方程、Schrödinger 方程等多变量问题常通过离散化转化为矩阵方程.

Example 二维 Poisson 方程

\[-\Delta u=f(x,y),\qquad f_{ij}=f(x_i,y_j).\]

Solution 在矩形网格上取步长 \(h_x,h_y\). 对内部点 \((x_i,y_j)\), 二阶中心差分给出

\[\left(\frac{2}{h_x^2}+\frac{2}{h_y^2}\right)U_{ij} -\frac{U_{i-1,j}+U_{i+1,j}}{h_x^2} -\frac{U_{i,j-1}+U_{i,j+1}}{h_y^2} =f_{ij}.\]

其中 \(U_{ij}\approx u(x_i,y_j)\) 替代了精确值. 已经舍去截断误差.

这就是五点差分格式: 每个内部点只与上下左右四个邻点耦合.

如果有混合导数, 将变为九点格式.

若进一步取 \(h_x=h_y=h\), 则

\[4U_{ij}-U_{i-1,j}-U_{i+1,j}-U_{i,j-1}-U_{i,j+1}=h^2 f_{ij}.\]

  • 内部格点可以投影为一维未知量. 对 \((m-1)\times(n-1)\) 个内部点取

    \[k=i+(m-1)(j-1)\]

    若把 \(h^2\) 乘到右端, 所有内部未知量组成向量 \(\mathbf{U}\), 离散方程写为

    \[\mathbf{A}\mathbf{U}=\mathbf{F}.\]

    五点差分格式给出分块三对角系数.

    例如对 \(4\times 4\) 二维网格的 \(3\times3\) 个内部点按行投影时,

    \[\begin{aligned} x_1&=U_{11},& x_2&=U_{21},& x_3&=U_{31},\\ x_4&=U_{12},& x_5&=U_{22},& x_6&=U_{32},\\ x_7&=U_{13},& x_8&=U_{23},& x_9&=U_{33}. \end{aligned}\]

    对这个例子

    \[\mathbf{A}= \begin{bmatrix} \mathbf{B}&-\mathbf{I}&0\\ -\mathbf{I}&\mathbf{B}&-\mathbf{I}\\ 0&-\mathbf{I}&\mathbf{B}\end{bmatrix}, \qquad \mathbf{B} = \begin{bmatrix} 4&-1&0\\ -1&4&-1\\ 0&-1&4 \end{bmatrix}\]

  • 边界条件进入矩阵或右端项.

    • Dirichlet 条件直接给出边界值;
    • Neumann 条件可用虚拟点或单边差分转化为内部点关系;
    • mixed 条件则把函数值和法向导数一起纳入离散方程.

    Example 左边界 \(x=x_0\) 的 Neumann 条件.

    Solution 这里定义 \(u_x(x_0,y)=g(x_0,y)\), 因而可用外部虚拟点表示为

    \[\frac{U_{1,j}-U_{-1,j}}{2h}=g(x_0,y_j) \quad\Longrightarrow\quad U_{-1,j}=U_{1,j}-2h\,g(x_0,y_j).\]

    边界点的五点方程

    \[4U_{0,j}-U_{-1,j}-U_{1,j}-U_{0,j+1}-U_{0,j-1} =h^2 f(x_0,y_j)\]

    消去虚拟点后得到

    \[4U_{0,j}-2U_{1,j}-U_{0,j+1}-U_{0,j-1} =h^2 f(x_0,y_j)-2h\,g(x_0,y_j).\]

    若给定的是外法向导数 \(\partial_nu=g_n\), 左边界上 \(g=-g_n\), 因而右端改为 \(h^2f+2hg_n\).

Example 求解

\[-\Delta u=f(x,y) =-2\cos \pi y+\pi^2(1+x)^2\cos \pi y,\qquad 0\lt x,y\lt 1,\]

边界为

\[\begin{aligned} u(x,0)&=(1+x)^2,& u(x,1)&=-(1+x)^2,\\ u(0,y)&=\cos\pi y,& u(1,y)&=4\cos\pi y, \end{aligned} \qquad h=\frac14.\]

精确解是 \(u(x,y)=(1+x)^2\cos\pi y\).

Solution 内部节点五点方程同前. 对于 Dirichlet 边界点

  • 在 \((x_1,y_1)=(1/4,1/4)\)

    \[4U_{1,1}-U_{0,1}-U_{2,1}-U_{1,0}-U_{1,2} =h^2 f\left(\frac14,\frac14\right),\]

    代入 \(U_{0,1}=\cos(\pi/4)\), \(U_{1,0}=(1+1/4)^2\) 得

    \[4U_{1,1}-U_{2,1}-U_{1,2} =h^2 f\left(\frac14,\frac14\right) +\cos\frac{\pi}{4}+\left(1+\frac14\right)^2.\]

  • \((x_2,y_1)=(1/2,1/4)\)

    \[-U_{1,1}+4U_{2,1}-U_{3,1}-U_{2,2} =h^2 f\left(\frac12,\frac14\right)+\left(1+\frac12\right)^2.\]

  • \((x_3,y_1)=(3/4,1/4)\)

    \[-U_{2,1}+4U_{3,1}-U_{3,2} =h^2 f\left(\frac34,\frac14\right) +4\cos\frac{\pi}{4}+\left(1+\frac34\right)^2.\]

其余六个节点也按相同五点格式列式, 代入已知边界值, 得到九元线性方程组. 因外力及边界关于 \(y=1/2\) 反对称, 离散解满足 \(U_{i,2}=0\)、\(U_{i,3}=-U_{i,1}\). 上面的前三式因此闭合为

\[\begin{pmatrix}4&-1&0\\-1&4&-1\\0&-1&4\end{pmatrix} \begin{pmatrix}U_{1,1}\\U_{2,1}\\U_{3,1}\end{pmatrix} =\begin{pmatrix}b_1\\b_2\\b_3\end{pmatrix},\]

其中 \(b_1,b_2,b_3\) 分别为前三式已给出的右端. 求解得全部内部值 (行对应 \(y_j\), 列对应 \(x_i\)):

\[(U_{i,j})_{j,i=1}^{3}\approx \begin{pmatrix} 1.1187788720&1.6123673475&2.1876760876\\ 0&0&0\\ -1.1187788720&-1.6123673475&-2.1876760876 \end{pmatrix}.\]

这是五点格式的离散解. 与精确连续解在九个内部节点处比较, 最大绝对误差约为 \(0.02216157\); 两者不应混为同一个解.

二维、三维问题的网格点数随维度快速增长, 很难单靠缩小格点间距提高精度.

实际计算要同时利用边界条件、对称性和物理量子数改进微分算符.

10.2 Fourier Derivative Operators and FFT(Fast Fourier Transform)

Fourier 变换可以用来数值求解导数. 在周期性边界条件下, 可把导数计算从实空间转换到 Fourier 空间.

考虑一维等间距格点, 均匀分布在 \([-L/2,L/2)\) 上

\[x_i=-\frac{L}{2}+i\Delta x =\left(i-\frac{N}{2}\right)\Delta x,\qquad i=0,1,\cdots,N-1,\]

其中 \(L=N\Delta x\), 以下取偶数 \(N\). 周期边界条件为 \(f(x+L)=f(x)\).

  • FFT 中的未折叠频率索引为

    \[k_n=\frac{2\pi n}{L},\qquad n=0,1,\cdots,N-1.\]

    标准 FFT 排列把 \(n\ge N/2\) 解释为负频率:2

    \[k_n= \begin{cases} \dfrac{2\pi n}{L}, & n=0,1,\cdots,N/2-1,\\[4pt] \dfrac{2\pi(n-N)}{L}, & n=N/2,\cdots,N-1, \end{cases} \qquad \Delta k=\frac{2\pi}{L}.\]

    其中 Nyquist 模 \(n=N/2\) 的正负频率在网格上不可区分. 实数三角插值采用对称的余弦 Nyquist 模: 奇数阶求导时该模的乘子取 0, 偶数阶时保留 \((\mathrm{i}k_{N/2})^m\).

  • DFT(discrete Fourier transform) 为

    \[\hat f(k_n)=\sum_{i=0}^{N-1} \exp(-\mathrm{i}k_nx_i)f(x_i) \Longleftrightarrow f(x_i)=\frac{1}{N}\sum_{n=0}^{N-1} \exp(\mathrm{i}k_nx_i)\hat f(k_n).\]

    连续 Fourier 变换中的求导对应乘以 \((\mathrm{i}k)^m\). 对有限格点, 以下等式给出由样本构造的三角插值函数 \(f_N\) 的导数, 并按上述规则处理 Nyquist 模; 只有当原函数属于该插值空间时, 它才等于原函数的精确导数. 一般光滑周期函数须以增加 \(N\) 检查谱近似的收敛. 因此

    \[f_N^{(m)}(x_i) =\frac{1}{N}\sum_{n=0}^{N-1} \exp(\mathrm{i}k_nx_i)\mu_n^{(m)} \sum_{j=0}^{N-1}\exp(-\mathrm{i}k_nx_j)f(x_j) =\sum_{j=0}^{N-1}D_{ij}^{(m)}f(x_j)\]

    其中

    \[\mu_n^{(m)}=\begin{cases} 0, & n=N/2,\quad m\text{ 为奇数},\\ (\mathrm{i}k_n)^m, & \text{其余情形}, \end{cases}\]

    \[D_{ij}^{(m)} =\frac{1}{N}\sum_{n=0}^{N-1} \exp[\mathrm{i}k_n(x_i-x_j)]\mu_n^{(m)}.\]

    • 一阶微分矩阵元为

      \[D_{ij}^{(1)} =\frac{1}{N}\sum_{n=0}^{N-1} \exp[\mathrm{i}k_n(x_i-x_j)]\mu_n^{(1)}.\]

      对上述实数三角插值, 将 Nyquist 乘子设为 0 后, 一阶微分矩阵为

      \[D_{ij}^{(1)} =-\frac{2\Delta k}{N} \sum_{n=1}^{N/2-1}n\sin[k_n\Delta x(i-j)].\]

    • 同理二阶微分矩阵元为

      \[D_{ij}^{(2)} =-\frac{2(\Delta k)^2}{N} \sum_{n=1}^{N/2-1}n^2\cos[k_n\Delta x(i-j)] -\frac{(\Delta k)^2}{N} \left(\frac{N}{2}\right)^2 \cos\!\left[\frac{N\Delta k\,\Delta x}{2}(i-j)\right].\]

Remark 对于 Nyquist 样本 \(f_i=(-1)^i\), 有 \(D^{(1)}\mathbf{f}=0\), 而 \(D^{(2)}\mathbf{f}=-(\pi/\Delta x)^2\mathbf{f}\). 因而此约定下 \(D^{(2)}\ne(D^{(1)})^2\), 二阶导数应直接乘以 \(-k^2\) 计算.3

Example 利用 FFT 求解波动型方程

\[\partial_t^2 u = C(x,t)\,\partial_x^2 u\]

Solution 时间上差分推进, 空间上用 FFT 做谱微分.

记 \(u_j^n \approx u(x_j,t_n)\), 每个时刻的流程为:

  1. 对 \(u^n(x)\) 作 FFT, 得到 \(\widehat u_k^n\);
  2. 在 Fourier 空间中乘以 \(-k^2\): \[\widehat{u_{xx}}(k,t_n) = -k^2\,\hat u_k^n;\]
  3. 逆 FFT 回到实空间, 得到 \(u_{xx}(x_j,t_n)\);
  4. 在实空间乘变系数: \[F_j^n = C(x_j,t_n)\,u_{xx}(x_j,t_n).\]

时间推进分两步:

  1. 构造第一层 \(u_j^1\): 由 Taylor 展开

    \[u_j^1 = u_j^0 + \Delta t\,v(x_j) + \frac{\Delta t^2}{2}F_j^0, \qquad v(x)=u_t(x,0);\]

  2. 对 \(n\ge 1\) 用中心差分

    \[u_j^{n+1} = 2u_j^n - u_j^{n-1} + \Delta t^2 F_j^n.\]

Proposition 对于常系数 \(C=c^2>0\), 每个 Fourier 模是频率 \(c|k|\) 的振子; 非零模保持有界的稳健步长条件为

\[c\Delta t\max_n|k_n|\lt2.\]

偶数网格下 \(\max|k_n|=\pi/\Delta x\), 因而

\[c\Delta t/\Delta x\lt2/\pi.\]

等号处最高频模的特征方程为 \((G+1)^2=0\); 除非该模的起步满足 \(u^1=-u^0\), 否则会出现 \((-1)^n(A+Bn)\) 的线性增长, 实际步长须留裕量.

Remark 变系数时须检查离散空间算符的谱, 不能直接把有限差分波动方程的 \(c\Delta t/\Delta x\le1\) 套到谱微分上.

因为 \(C(x,t)\) 依赖空间位置, 方程在 Fourier 空间中不能完全逐模对角化, 因此采用"频域求导, 实空间乘系数"的混合实现. FFT 最适合周期边界条件; 非周期问题须选择匹配边界条件的基, 例如齐次 Dirichlet/Neumann 条件可分别用正弦/余弦展开.

10.3 Basis Expansion Methods *

基矢展开法的计算目标是把微分方程或本征值问题转化为有限维矩阵问题.

以 Schrödinger 方程 \(\hat H|\psi\rangle=E|\psi\rangle\) 为例. 选一组正交归一基 \(\{|\phi_j\rangle\}\), 展开

\[|\psi\rangle=\sum_j c_j|\phi_j\rangle.\]

代入本征值方程并左乘 \(\langle\phi_i|\), 得

\[\sum_j H_{ij}c_j=Ec_i, \qquad H_{ij}=\langle\phi_i|\hat H|\phi_j\rangle.\]

若基不正交归一, 则须保留重叠矩阵 \(S_{ij}=\langle\phi_i|\phi_j\rangle\), 求解 \(\sum_jH_{ij}c_j=E\sum_jS_{ij}c_j\), 并按 \(\mathbf{c}^\dagger\mathbf{S}\mathbf{c}=1\) 归一化.

计算流程:

  1. 选择基函数 \(\{|\phi_j\rangle\}\), 例如谐振子基、平面波基、DVR(discrete variable representation) 基等.
  2. 计算矩阵元 \(H_{ij}\).
  3. 截断到有限维矩阵.
  4. 对角化矩阵, 得到近似本征值和展开系数.

10.4 DVR(Discrete Variable Representation) *

DVR(discrete variable representation) 构造一组特殊正交基, 使波函数展开系数正比于格点值, 并把内积写成加权求和. 设投影算符为 \(P\), 网格点为 \(x_\alpha\), 定义未归一化 DVR 态

\[|\Delta_\alpha\rangle=P|x_\alpha\rangle.\]

若

\[\langle\Delta_\alpha|\Delta_\beta\rangle=N_\alpha\delta_{\alpha\beta},\]

则归一化 DVR 基为

\[|F_\alpha\rangle=\frac{1}{\sqrt{N_\alpha}}|\Delta_\alpha\rangle.\]

核心性质:

\[\langle F_\alpha|F_\beta\rangle=\delta_{\alpha\beta}, \qquad F_\beta(x_\alpha)=\sqrt{N_\alpha}\delta_{\alpha\beta}.\]

以下假设 \(P\) 是正交投影, \(\{|F_\alpha\rangle\}\) 张成其像空间, 且 \(P|\psi\rangle=|\psi\rangle\) (内积中的 \(\phi\) 也属于该子空间). 因此

\[\langle F_\alpha|\psi\rangle=\frac{\psi(x_\alpha)}{\sqrt{N_\alpha}}, \qquad |\psi\rangle=\sum_\alpha\frac{\psi(x_\alpha)}{\sqrt{N_\alpha}}|F_\alpha\rangle.\]

内积可写为

\[\langle\phi|\psi\rangle =\sum_\alpha\frac{1}{N_\alpha}\phi^*(x_\alpha)\psi(x_\alpha) =\sum_\alpha\omega_\alpha\phi^*(x_\alpha)\psi(x_\alpha), \qquad \omega_\alpha=N_\alpha^{-1}.\]

从正交基构造 DVR: 若

\[P=\sum_n|\phi_n\rangle\langle\phi_n|, \qquad |\Delta_\alpha\rangle=P|x_\alpha\rangle =\sum_n|\phi_n\rangle\phi_n^*(x_\alpha),\]

则

\[\Delta_\alpha(x)=\sum_n\phi_n(x)\phi_n^*(x_\alpha), \qquad \langle\Delta_\alpha|\Delta_\beta\rangle =\sum_n\phi_n(x_\alpha)\phi_n^*(x_\beta).\]

若取 \(\phi_n(x)=\sqrt{W(x)}P_n(x)\), 其中 \(P_n\) 在权重 \(W\) 下正交归一, 保留 \(n=0,\cdots,M-1\), 且 \(x_\alpha\) 取 \(P_M\) 的 \(M\) 个零点. 记对应的 Gauss 权为 \(\lambda_\alpha\), 即

\[\int W(x)Q(x)\,\mathrm{d}x=\sum_\alpha\lambda_\alpha Q(x_\alpha),\qquad \deg Q\le2M-1.\]

由于 \(\sum_\alpha\lambda_\alpha P_n(x_\alpha)P_m(x_\alpha)=\delta_{nm}\), 方阵的正交性给出 \(\sum_{n=0}^{M-1}P_n(x_\alpha)P_n(x_\beta)=\delta_{\alpha\beta}/\lambda_\alpha\). 因而 \(N_\alpha=W(x_\alpha)/\lambda_\alpha\), 位置内积的求积权为 \(w_\alpha=\omega_\alpha=\lambda_\alpha/W(x_\alpha)\), 并有

\[\sum_n\phi_n(x_\alpha)\phi_n^*(x_\beta)=N_\alpha\delta_{\alpha\beta}.\]

归一化 DVR 基可写为

\[F_\alpha(x)=\sqrt{w_\alpha}\sum_n\sqrt{W(x)W(x_\alpha)}P_n(x)P_n(x_\alpha).\]

在 DVR 基下, 势能矩阵近似对角:

\[\langle F_\alpha|V|F_\beta\rangle \approx\sum_kw_kF_\alpha^*(x_k)V(x_k)F_\beta(x_k) =\delta_{\alpha\beta}V(x_\alpha).\]

平面波 DVR 例子: 取动量截断 \(p\in[-p_0,p_0]\),

\[P=\int_{-p_0}^{p_0}|p\rangle\langle p|\,\mathrm{d}p, \qquad x_\alpha=\frac{\pi\hbar}{p_0}\alpha,\qquad \alpha\in\mathbb{Z}.\]

未归一化基为

\[|\Delta_\alpha\rangle =\int_{-p_0}^{p_0}\mathrm{d}p\,|p\rangle \frac{\exp(-\mathrm{i}px_\alpha/\hbar)}{\sqrt{2\pi\hbar}}.\]

内积为

\[\langle\Delta_\alpha|\Delta_\beta\rangle =\int_{-p_0}^{p_0}\mathrm{d}p\, \frac{\exp[\mathrm{i}p(x_\alpha-x_\beta)/\hbar]}{2\pi\hbar} =\frac{p_0}{\pi\hbar} \frac{\sin[p_0(x_\alpha-x_\beta)/\hbar]}{p_0(x_\alpha-x_\beta)/\hbar} =\frac{p_0}{\pi\hbar}\delta_{\alpha\beta}.\]

归一化基为

\[|F_\alpha\rangle=\sqrt{\frac{\pi\hbar}{p_0}}|\Delta_\alpha\rangle.\]

位置表象:

\[F_\alpha(x)=\frac{1}{\sqrt{\pi\hbar p_0}} \frac{\sin[p_0(x-x_\alpha)/\hbar]}{(x-x_\alpha)/\hbar} =\sqrt{\frac{p_0}{\pi\hbar}}\, \mathrm{sinc}\left(\frac{p_0(x-x_\alpha)}{\hbar}\right).\]

这里 \(\mathrm{sinc}(z)\equiv\sin z/z\), 并以连续极限定义 \(\mathrm{sinc}(0)=1\). 以下展开针对动量支撑在 \([-p_0,p_0]\) 内的带限波函数; 无限网格的正交与完备性不等于有限截断后可精确表示任意波函数.

动量表象与波函数展开:

\[\langle p|F_\alpha\rangle=\frac{1}{\sqrt{2p_0}}\exp(-\mathrm{i}px_\alpha/\hbar), \qquad |\psi\rangle=\sum_\alpha |F_\alpha\rangle\sqrt{\frac{\pi\hbar}{p_0}}\psi(x_\alpha).\]

10.5 FEM(Finite-Element Method) *

有限元法适合复杂几何边界. 基本流程:

  1. 将区域剖分为互不重叠的小单元.
  2. 在每个单元上用低阶多项式近似待求场量.
  3. 用变分原理或加权残量法把微分方程离散为代数方程组.
  4. 组装全局矩阵并求解.

Delaunay 三角剖分常用于生成三角网格: 最大化最小角、避免瘦长三角形、三角形并集为点集凸包.

10.5.1 Variational Formulation

考虑二维 Poisson 方程

\[\nabla^2\phi(\mathbf{r})=f(\mathbf{r}), \qquad \mathbf{r}\in\Omega.\]

本节沿用 \(\nabla^2\phi=f\) 的符号约定. 以下先取给定 Dirichlet 边界, 使允许变分在边界上为零. 对应变分泛函可写为

\[J[\phi]=\int_\Omega\left[\frac12(\nabla\phi)^2+f\phi\right]\,\mathrm{d}^2r, \qquad \delta J=0.\]

在三角单元上, 线性插值为

\[\phi(x,y)=a_i+b_ix+c_iy,\]

二次插值可写为

\[\phi(x,y)=a_i+b_ix+c_iy+d_ixy+e_ix^2+f_iy^2.\]

10.5.2 Area Coordinates

以下将顶点 \(a,b,c\) 按逆时针顺序编号, 并令 \(A>0\) 为三角形面积, 从而后面的行列式为 \(2A\). 自然坐标 \(\xi_a,\xi_b,\xi_c\) 定义为面积比:

\[\xi_a=\frac{A_a}{A}, \qquad \xi_b=\frac{A_b}{A}, \qquad \xi_c=\frac{A_c}{A}, \qquad \xi_a+\xi_b+\xi_c=1.\]

顶点处满足 \(\xi_a=1,\xi_b=\xi_c=0\) 及循环关系.

Finite-element natural coordinates in a triangular element

坐标变换:

\[\begin{pmatrix}1\\x\\y\end{pmatrix} = \begin{pmatrix}1&1&1\\x_a&x_b&x_c\\y_a&y_b&y_c\end{pmatrix} \begin{pmatrix}\xi_a\\\xi_b\\\xi_c\end{pmatrix}.\]

逆变换:

\[\begin{pmatrix}\xi_a\\\xi_b\\\xi_c\end{pmatrix} =\frac{1}{2A} \begin{pmatrix} x_by_c-x_cy_b&y_b-y_c&x_c-x_b\\ x_cy_a-x_ay_c&y_c-y_a&x_a-x_c\\ x_ay_b-x_by_a&y_a-y_b&x_b-x_a \end{pmatrix} \begin{pmatrix}1\\x\\y\end{pmatrix},\]

其中

\[2A= \begin{vmatrix}1&1&1\\x_a&x_b&x_c\\y_a&y_b&y_c\end{vmatrix}.\]

单元内函数值:

\[\phi(x,y)=\xi_a\phi_a+\xi_b\phi_b+\xi_c\phi_c.\]

面积坐标积分公式:

\[\int_A\xi_a^k\xi_b^l\xi_c^m\,\mathrm{d}\Omega =2A\frac{k!l!m!}{(k+l+m+2)!}.\]

10.5.3 Element Matrices and Assembly

代入线性插值得

\[\nabla\phi=\nabla\xi_a\phi_a+\nabla\xi_b\phi_b+\nabla\xi_c\phi_c.\]

于是单元泛函化为

\[J[\phi] =\frac12\boldsymbol{\phi}^{\mathsf{T}}\mathbf{K}\boldsymbol{\phi} +\mathbf{r}^{\mathsf{T}}\boldsymbol{\phi},\]

其中

\[\boldsymbol{\phi}=(\phi_a,\phi_b,\phi_c)^{\mathsf{T}}, \qquad K_{ij}=\int_A\nabla\xi_i\cdot\nabla\xi_j\,\mathrm{d}\Omega, \qquad r_i=\int_A f\xi_i\,\mathrm{d}\Omega.\]

变分后得到线性方程组

\[\mathbf{K}\boldsymbol{\phi}=-\mathbf{r}.\]

因为本节采用 \(\nabla^2\phi=f\), 泛函中源项的符号为 \(+f\phi\), 驻值条件给出 \(\mathbf{K}\boldsymbol{\phi}+\mathbf{r}=0\). 若改写为 \(-\nabla^2\phi=f\), 则泛函应改用 \(-f\phi\), 才得到 \(\mathbf{K}\boldsymbol{\phi}=\mathbf{r}\). 全局组装后还须施加边界条件; 单个三角元的刚度矩阵有常数零模, 不能单独求逆.

自然坐标梯度:

\[\nabla\xi_a=\frac{1}{2A}(y_b-y_c,x_c-x_b),\]\[\nabla\xi_b=\frac{1}{2A}(y_c-y_a,x_a-x_c), \qquad \nabla\xi_c=\frac{1}{2A}(y_a-y_b,x_b-x_a).\]

令

\[\mathbf{b}=(y_b-y_c,y_c-y_a,y_a-y_b)^{\mathsf{T}}, \qquad \mathbf{c}=(x_c-x_b,x_a-x_c,x_b-x_a)^{\mathsf{T}}.\]

单元刚度矩阵为

\[\mathbf{K}=\frac{1}{4A}(\mathbf{b}\mathbf{b}^{\mathsf{T}}+\mathbf{c}\mathbf{c}^{\mathsf{T}}).\]

等价地, 记 \(y_{ab}=y_a-y_b\), \(y_{bc}=y_b-y_c\), \(y_{ca}=y_c-y_a\), 以及 \(x_{ab}=x_a-x_b\), \(x_{bc}=x_b-x_c\), \(x_{ca}=x_c-x_a\), 则

\[\mathbf{K}=\frac{1}{4A} \begin{pmatrix} y_{bc}^2+x_{bc}^2&y_{bc}y_{ca}+x_{bc}x_{ca}&y_{bc}y_{ab}+x_{bc}x_{ab}\\ y_{ca}y_{bc}+x_{ca}x_{bc}&y_{ca}^2+x_{ca}^2&y_{ca}y_{ab}+x_{ca}x_{ab}\\ y_{ab}y_{bc}+x_{ab}x_{bc}&y_{ab}y_{ca}+x_{ab}x_{ca}&y_{ab}^2+x_{ab}^2 \end{pmatrix}.\]

若源项也用线性插值

\[f(\mathbf{r})=\xi_af_a+\xi_bf_b+\xi_cf_c,\]

则载荷向量为

\[r_a=2A\left(\frac{f_a}{12}+\frac{f_b}{24}+\frac{f_c}{24}\right),\]\[r_b=2A\left(\frac{f_b}{12}+\frac{f_c}{24}+\frac{f_a}{24}\right), \qquad r_c=2A\left(\frac{f_c}{12}+\frac{f_a}{24}+\frac{f_b}{24}\right).\]

对每个三角元计算 \(\mathbf{K}\) 和 \(\mathbf{r}\), 再组装成全局刚度矩阵和载荷向量, 最后求解全局线性方程组.

11 Dynamical Equations

11.1 One-Dimensional Diffusion Equation

一维扩散方程为

\[\frac{\partial u}{\partial t} =\frac{\partial^2u}{\partial x^2}, \qquad 0\lt x\lt L.\]

初始条件为 \(u(x,0)=g(x)\). 有限区间还须给定边界条件; 以下差分矩阵取齐次 Dirichlet 边界 \(u(0,t)=u(L,t)=0\). 为说明扩散核, 先单独考虑全空间 \(x\in\mathbb{R}\): 单位点源 \(\delta\) 函数的基本解为

\[w(x,t)=\frac{1}{\sqrt{4\pi t}}\exp\left(-\frac{x^2}{4t}\right).\]

在全空间且初值满足积分收敛条件时, 解为基本解与初始条件的卷积

\[u(x,t)=\int_{-\infty}^{\infty}w(x-x',t)g(x')\,\mathrm{d}x'.\]

有限区间的齐次 Dirichlet 问题则用满足边界条件的正弦展开:

\[u(x,t)=\sum_{n=1}^{\infty}b_n\sin\frac{n\pi x}{L}\,\mathrm{e}^{-(n\pi/L)^2t},\qquad b_n=\frac{2}{L}\int_0^L g(x)\sin\frac{n\pi x}{L}\,\mathrm{d}x.\]

例如 \(g\in L^2(0,L)\) 时, 正弦基的完备性保证初值在 \(L^2\) 意义恢复, 而 \(t>0\) 的指数衰减保证解平滑. 下文回到该有限区间问题.

11.1.1 Finite-Difference Method

记 \(u_{i,j}=u(x_i,t_j)\), 时间差分和空间差分分别为

\[\frac{\partial u}{\partial t} \approx\frac{u_{i,j+1}-u_{i,j}}{\Delta t}, \qquad \frac{\partial^2u}{\partial x^2} \approx\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{\Delta x^2}.\]

定义 \(\alpha=\frac{\Delta t}{\Delta x^2}\), 得到递推公式

\[u_{i,j+1}=\alpha u_{i-1,j}+(1-2\alpha)u_{i,j}+\alpha u_{i+1,j}.\]
  • 对 Fourier 模 \(u_{i,j}=G^j\mathrm{e}^{\mathrm{i}i\vartheta}\), 放大因子为 \(G=1-4\alpha\sin^2(\vartheta/2)\). 要求所有模式满足 \(|G|\le1\), 即 \(0\le\alpha\le1/2\). 若扩散系数为 \(D\), 应改用 \(\alpha=D\Delta t/\Delta x^2\).
  • 第一步为 \(u_{i,1}=\alpha g(x_{i-1})+(1-2\alpha)g(x_i)+\alpha g(x_{i+1})\).

11.1.2 Crank-Nicolson Scheme

Crank-Nicolson 格式中, 时间导数采用向后差分, 空间二阶导数取两个时间层的加权平均:

\[\frac{\partial u}{\partial t} \approx\frac{u_{i,j}-u_{i,j-1}}{\Delta t}, \qquad \frac{\partial^2u}{\partial x^2} \approx \theta\left.\frac{\partial^2u}{\partial x^2}\right|_{t=t_j} +(1-\theta)\left.\frac{\partial^2u}{\partial x^2}\right|_{t=t_{j-1}}.\]

离散后得到

\[\theta\frac{u_{i-1,j}-2u_{i,j}+u_{i+1,j}}{\Delta x^2} +(1-\theta)\frac{u_{i-1,j-1}-2u_{i,j-1}+u_{i+1,j-1}}{\Delta x^2} =\frac{u_{i,j}-u_{i,j-1}}{\Delta t}.\]

\(\theta=0\) 时为显式格式; \(\theta=1\) 时为全隐式格式; \(\theta=1/2\) 时为 Crank-Nicolson 格式.

取 \(\theta=1/2\), 并令 \(\alpha=\Delta t/\Delta x^2\), 整理得

\[-\alpha u_{i-1,j}+2(1+\alpha)u_{i,j}-\alpha u_{i+1,j} =\alpha u_{i-1,j-1}+2(1-\alpha)u_{i,j-1}+\alpha u_{i+1,j-1}.\]

可以记为

\[(2\mathbf{I}+\alpha\mathbf{B})\mathbf{u}_j =(2\mathbf{I}-\alpha\mathbf{B})\mathbf{u}_{j-1}, \qquad \mathbf{B}= \begin{pmatrix} 2&-1&0&\cdots&0\\ -1&2&-1&\cdots&0\\ 0&\ddots&\ddots&\ddots&\vdots\\ \vdots&\ddots&-1&2&-1\\ 0&\cdots&0&-1&2 \end{pmatrix},\]

则递推关系

\[\mathbf{u}_j =(2\mathbf{I}+\alpha\mathbf{B})^{-1}(2\mathbf{I}-\alpha\mathbf{B})\mathbf{u}_{j-1}.\]

实现时求解三对角方程组, 不显式构造逆矩阵. 对 Fourier 模, Crank-Nicolson 放大因子为 \(G=[1-2\alpha\sin^2(\vartheta/2)]/[1+2\alpha\sin^2(\vartheta/2)]\), 任意 \(\alpha\ge0\) 均满足 \(|G|\le1\). 无条件稳定不等于任意步长都准确: 大步长使高频模式趋近 \(G=-1\), 可能出现缓慢衰减的数值振荡. 光滑解的全局误差为 \(O(\Delta t^2+\Delta x^2)\).

11.2 Wave Equation

11.2.1 One-Dimensional Wave Equation

一维波动方程为

\[\frac{\partial^2u}{\partial x^2} =\frac{\partial^2u}{\partial t^2}.\]

定解条件为

\[u(x,0)=g(x), \qquad \left.\frac{\partial u}{\partial t}\right|_{t=0}=0, \qquad u(0,t)=u(1,t)=0.\]

  • 中心差分为

    \[\frac{\partial^2u}{\partial x^2} \approx\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{\Delta x^2}, \qquad \frac{\partial^2u}{\partial t^2} \approx\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{\Delta t^2}.\]

    得到递推公式

    \[u_{i,j+1}=2u_{i,j}-u_{i,j-1} +\frac{\Delta t^2}{\Delta x^2}(u_{i+1,j}-2u_{i,j}+u_{i-1,j}).\]

  • 初速度条件给出自然边界条件

    \[\frac{u_{i,1}-u_{i,-1}}{2\Delta t}=0, \qquad u_{i,1}=u_{i,-1}.\]

    代入递推公式可得

    \[u_{i,1}=g(x_i) +\frac{\Delta t^2}{2\Delta x^2} \left[g(x_{i+1})-2g(x_i)+g(x_{i-1})\right].\]

一维格式的 Fourier 放大因子满足 \(G^2-2[1-2r^2\sin^2(\vartheta/2)]G+1=0\), 其中 \(r=c\Delta t/\Delta x\). 常系数问题的稳定步长要求 \(r\le1\), 本文取 \(c=1\). 该约束是所有网格模式的要求, 不能仅靠某个平滑初值下的短时间测试判断稳定性.

11.2.2 Two-Dimensional Wave Equation

二维波动方程为

\[\frac{\partial^2u}{\partial x^2} +\frac{\partial^2u}{\partial y^2} =\frac{\partial^2u}{\partial t^2}, \qquad (x,y)\in(0,1)\times(0,1).\]

定解条件为初值 \(u(x,y,0)=g(x,y)\), 初速度 \(u_t(x,y,0)=0\), 边界取零值.

  • 离散化取

    \[t_l=l\Delta t, \qquad x_i=i\Delta x, \qquad y_j=j\Delta y, \qquad \Delta x=\Delta y=h.\]

    中心差分为

    \[u_{xx}\approx\frac{u_{i+1,j}^{l}-2u_{i,j}^{l}+u_{i-1,j}^{l}}{h^2}, \quad u_{yy}\approx\frac{u_{i,j+1}^{l}-2u_{i,j}^{l}+u_{i,j-1}^{l}}{h^2}, \quad u_{tt}\approx \frac{u_{i,j}^{l+1}-2u_{i,j}^{l}+u_{i,j}^{l-1}}{\Delta t^2}.\]

    代入得到递推公式

    \[u_{i,j}^{l+1}=2u_{i,j}^{l}-u_{i,j}^{l-1} +\frac{\Delta t^2}{h^2} \left(u_{i+1,j}^{l}+u_{i-1,j}^{l}+u_{i,j+1}^{l}+u_{i,j-1}^{l}-4u_{i,j}^{l}\right).\]

  • 零初速给出 \(u_{i,j}^{-1}=u_{i,j}^{1}\), 因而

    \[u_{i,j}^{1}=g(x_i,y_j) +\frac{\Delta t^2}{2h^2} \left[g(x_{i+1},y_j)+g(x_{i-1},y_j)+g(x_i,y_{j+1})+g(x_i,y_{j-1})-4g(x_i,y_j)\right].\]

二维格式同理要求 \(c^2\Delta t^2(\Delta x^{-2}+\Delta y^{-2})\le1\). 本节 \(c=1\), 等间距网格下即 \(\Delta t/h\le1/\sqrt2\).

11.3 Time Evolution in Quantum Mechanics *

取 \(\hbar=m=1\) 的单位, 考虑 Schrödinger 方程

\[\mathrm{i}\frac{\partial}{\partial t}\psi(x,t) =\left[-\frac{1}{2}\frac{\partial^2}{\partial x^2}+V(x)+\epsilon(t)x\right]\psi(x,t).\]

空间二阶导数用中心差分:

\[\frac{\partial^2\psi_j}{\partial x^2} \approx \frac{\psi_{j+1}-2\psi_j+\psi_{j-1}}{\Delta x^2}.\]

于是离散 Hamiltonian 为三对角矩阵:

\[(\mathbf{H}\boldsymbol{\psi})_j =-\frac{1}{2\Delta x^2} (\psi_{j+1}-2\psi_j+\psi_{j-1}) +\left[V(x_j)+\epsilon(t)x_j\right]\psi_j.\]

即

\[H_{jj}=\frac{1}{\Delta x^2}+V(x_j)+\epsilon(t)x_j, \qquad H_{j,j\pm1}=-\frac{1}{2\Delta x^2}.\]

时间推进可用 Crank-Nicolson 形式:

\[\left(\mathbf{I}+\frac{\mathrm{i}\Delta t}{2}\mathbf{H}^{n+1/2}\right) \boldsymbol{\psi}^{n+1} = \left(\mathbf{I}-\frac{\mathrm{i}\Delta t}{2}\mathbf{H}^{n+1/2}\right) \boldsymbol{\psi}^n.\]

因此每一步只需求解一个复数三对角线性方程组. 左右两侧矩阵的非零元为

\[U_{jj}^{\pm} =1\pm\frac{\mathrm{i}\Delta t}{2} \left[\frac{1}{\Delta x^2}+V(x_j)+\epsilon(t_{n+1/2})x_j\right], \qquad U_{j,j\pm1}^{\pm} =\mp\frac{\mathrm{i}\Delta t}{4\Delta x^2}.\]

其中左端使用 \(U^+\), 右端使用 \(U^-\).

当势能为实数, 且边界离散保持 \(\mathbf H^{n+1/2}\) 为 Hermitian 矩阵时, 该推进保持离散范数. 令 \(a=\Delta t/2\), 每个实特征值 \(E\) 对应的推进因子为 \((1-\mathrm{i}aE)/(1+\mathrm{i}aE)\), 模长为 1; 正交特征向量展开因而给出 \(\Delta x\sum_j|\psi_j^{n+1}|^2=\Delta x\sum_j|\psi_j^n|^2\). 此结论不消除时间离散的相位误差, 实际计算还应检查步长收敛及线性求解残差.

References

  1. L. F. Shampine and M. W. Reichelt, The MATLAB ODE Suite, §3, 固定步长 BDF 的阶数与稳定性. ↩︎

  2. NumPy, numpy.fft.fftfreq, 偶数长度的标准频率排列. ↩︎

  3. Steven G. Johnson, Notes on FFT-based differentiation, 三角插值与 Nyquist 模的奇偶阶求导; SciPy, fftpack.diff, 奇数阶 Nyquist 乘子约定. ↩︎