Computational Physics

Episode 1. 计算方法

Computational Physics

REF：计算方法/scherer computational physics/Hoffman 计算物理学/丁泽军讲义

涉及内容：经典力学问题/电磁学/QM本征值问题/统计力学相变/随机过程/二次量子化/严格对角化/开放系统耗散问题/固体物理/随机矩阵

Review：数值方法

插值

$n$ $(x_1,y_1),\cdots,(x_n,y_n)$ $f(x)$ $f(x_i)=y_i$ $e_n(x_i)=0$ ，显然误差满足

\begin{matrix} (1) & e_{n} (x) \propto \prod_{i}^{n} (x - x_{i}) \end{matrix}

$n$ $f(x)$ $n$ 的多项式，形式上

\begin{matrix} (2) & f (x) = a_{0} + a_{1} x + \dots + a_{n - 1} x^{n - 1} + e_{n} (x) \end{matrix}

直接的想法是求解方程组

\begin{matrix} (3) & y_{i} = f (x_{i}) = \sum_{j = 0}^{n - 1} a_{j} x_{i}^{j} for all i = 1, \dots, n \end{matrix}

Lagrange插值

插值的用意：

计算积分

\begin{matrix} (4) & I = \int_{a}^{b} f (x) d x \end{matrix}

$I=h\sum_iy_i$ 。外推这个表达式，有一般形式

\begin{matrix} (5) & I = (b - a) \sum_{i} y_{i} w_{i} \end{matrix}

$w_i$ $\sum_iw_i=1$ $\sim\mathcal O(h)$ $\sim\mathcal O(h^2)$

\begin{matrix} (6) & I_{i} = [f (x_{i} + h) - f (x_{i})] h / 2 ≃ f^{'} (x_{i}) h^{2} / 2 \end{matrix}

实际上对应Taylor展开得不同阶项

\begin{matrix} (7) & f (x_{i} + z) = f (x_{i}) + f^{'} (x_{i}) z + f^{″} (x_{i}) z^{2} / 2! + \dots \end{matrix}

$f(x)$ $\sim\mathcal O(h^3)$ $\{y_i(x_i):i=1,\cdots,n\}$

\begin{matrix} (8) & I = (x_{2} - x_{1}) (f_{1} + f_{2}) / 2 + (x_{3} - x_{2}) (f_{2} + f_{3}) / 2 + (x_{4} - x_{3}) (f_{3} + f_{4}) / 2 + \dots \end{matrix}

数据点等距的情况下

\begin{matrix} (9) & I = \frac{h}{2} (f_{1} + 2 f_{2} + \dots + 2 f_{n - 1} + f_{n}) \end{matrix}

这种方法实际上用两点插值

\begin{matrix} (10) & f (x) = f (a) + \frac{f (b) - f (a)}{b - a} (x - a) \end{matrix}

$f(x)$ 完成了解析积分运算。

进一步，考虑Lagrange插值

\begin{matrix} (11) & f (x) = \sum_{i} w_{i} l_{i} (x) = \sum_{i} y_{i} \frac{\prod_{j \neq i} (x - x_{j})}{\prod_{j \neq i} (x_{i} - x_{j})} \end{matrix}

考虑等距数据点，积分表示为

\begin{matrix} (12) & I = \sum_{i}^{n} \int_{x_{1}}^{x_{1} + (n - 1) h} y_{i} l_{i} (x) d x = \sum_{i}^{n} y_{i} \int_{0}^{n - 1} h \frac{\prod_{j \neq i} (x^{'} - x_{j}^{'})}{\prod_{j \neq i} (x_{i}^{'} - x_{j}^{'})} d x^{'} = h \sum_{i}^{n} y_{i} w_{i}^{(n)} \end{matrix}

$x=hx'$ ，最后的权重定义为

\begin{matrix} (13) & w_{i}^{(n)} := \int_{0}^{n - 1} \frac{\prod_{j \neq i} (x^{'} - x_{j}^{'})}{\prod_{j \neq i} (x_{i}^{'} - x_{j}^{'})} d x^{'} \end{matrix}

可直接查表。

注意到Lagrange函数的性质

\begin{matrix} (14) & l_{i} (x_{j}) = δ_{i j} \end{matrix}

定义

\begin{matrix} (15) & l (x) = \sum_{i = 1}^{n} l_{i} (x) - 1 \end{matrix}

可以验证

\begin{matrix} (16) & \sum_{i} l_{i} (x) = 1 \end{matrix}

同时也满足权重函数的归一性。

Newton-Cotes插值

对于Lagrange插值，增加一个数据点时，插值函数要完全重新计算；相对的Newton-Cotes插值在处理上简单很多，只需在插值多项式上额外增加一项即可。

考虑如下形式的插值多项式

\begin{matrix} (17) & f (x) = b_{0} + b_{1} (x - x_{1}) + b_{2} (x - x_{1}) (x - x_{2}) + \dots + b_{n} (x - x_{n}) \dots (x - x_{n}) \end{matrix}

$x=x_i$ ，分别得到

\begin{matrix} (18) & \begin{array}{r} b_{0} = f (x_{1}), b_{1} = \frac{f (x_{1}) - f (x_{2})}{x_{1} - x_{2}} := f [x_{1}, x_{2}], b_{2} = \frac{f [x_{1}, x_{2}] - f [x_{2}, x_{3}]}{x_{1} - x_{3}} := f [x_{1}, x_{2}, x_{3}], \dots \end{array} \end{matrix}

其他插值方法

切比雪夫插值，考虑一个例子

\begin{matrix} (19) & \int_{- 1}^{1} f (t) d t = \int_{0}^{π} f (\cos t) \sin t d t \end{matrix}

其中

\begin{matrix} (20) & f (\cos t) = \frac{1}{2 M} C_{0} + \frac{1}{M} \sum_{j} C_{j} \cos (j t) + \frac{1}{2 M} C_{M} \cos (M t) \end{matrix}

$j=0,\cdots,M$ ，系数可以通过FFT计算。

Runge-Kutta方法

从一个ODE出发

\begin{matrix} (21) & d y / d x = f (x, y) \end{matrix}

$y_{n+1}-y_n=\int_{x_n}^{x_{n+!}}f(x,y)\dd{x}$ 。一个trivial的处理方式是

\begin{matrix} (22) & \begin{aligned} y_{n + 1} & = y_{n} + \int_{x_{n}}^{x_{n + 1}} f (x, y) d x \\ ≃ y_{n} + (x_{n + 1} - x_{n}) f (x_{n}, y_{n}) + O (h^{2}) \\ ≃ y_{n} + \frac{x_{n + 1} - x_{n}}{2} [f (x_{n}, y_{n}) + f (x_{n + 1}, y_{n + 1})] + O (h^{3}) \end{aligned} \end{matrix}

$y_{n+1}$ $y_{n+1}$ $\tilde y_{n+1}$ $y_{n+1}$ 的值。

换个角度考虑上式的积分

\begin{matrix} (23) & \begin{aligned} y_{n + 1} & = y_{n} + \int_{0}^{h} f (x_{n} + z, y (x_{n} + z)) d z \\ = y_{n} + \int_{0}^{h} f (x_{n} + z, y_{n} + y_{n}^{'} z + y_{n}^{″} z^{2} / 2) d z \\ = y_{n} + \int_{0}^{h} f + f_{x} z + f_{x x} z^{2} / 2 + f_{y} (y_{n}^{'} z + y_{n}^{″} z^{2} / 2) + \dots d z \\ = y_{n} + f (x_{n}, y_{n}) h + f_{x} h^{2} / 2 + f_{y} y_{n}^{'} h^{2} / 2 + O (h^{3}) \\ = y_{n} + f (x_{n}, y_{n}) h + h^{2} [f_{x} (x_{n}, y_{n}) + f_{y} (x_{n}, y_{n}) f (x_{n}, y_{n})] / 2 + O (h^{3}) \end{aligned} \end{matrix}

$f(x,y)$ 的各阶偏导数，实际仍然不便利。如何避开导数的运算？Runge-Kutta方法在保证精度的情况下，避开了导数的运算（证明较复杂）。

Runge-Kutta 23

由于

\begin{matrix} (24) & \begin{aligned} y_{n + 1} & = y_{n} + f (x_{n}, y_{n}) h + h^{2} [f_{x} (x_{n}, y_{n}) + f_{y} (x_{n}, y_{n}) f (x_{n}, y_{n})] / 2 + O (h^{3}) \\ = y_{n} \end{aligned} \end{matrix}

则

\begin{matrix} (25) & y_{n + 1} = y_{n} + h [a f (x_{n}, y_{n}) + b f (x_{n} + h, y_{n} + h)] + O (h^{3}) \end{matrix}

这个计算方法有3阶精度。

Runge-Kutta 45

令

\begin{matrix} (26) & \begin{aligned} k_{1} & = f (x_{n}, y_{n}) \\ k_{2} & = f (x_{n} + h / 2, y_{n} + h k_{1} / 2) \\ k_{3} & = f (x_{n} + h / 2, y_{n} + h k_{2} / 2) \\ k_{4} & = f (x_{n} + h, y_{n} + k 3) \end{aligned} \end{matrix}

则

\begin{matrix} (27) & y_{n + 1} = y_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}) + O (h^{5}) \end{matrix}

上述数值方法也称为ODE45。

求解二阶ODE

物理学问题中常常有二阶微分方程，例如单摆的运动方程

\begin{matrix} (28) & m l \ddot{θ} = m g \sin θ \end{matrix}

$v=\dot \theta$ ，reduce至一阶：

\begin{matrix} (29) & \begin{matrix} \frac{d}{d t} (\begin{matrix} θ \\ v \end{matrix}) = (\begin{matrix} v \\ g \sin θ / l \end{matrix}) \end{matrix} \end{matrix}

$y=\left( \begin{matrix} y_1\\y_2 \end{matrix} \right)$ 的向量情形。

HW ：Lagrange力学模拟问题

a.双摆问题

$U=-A/r^\alpha$

c. 竖杆倾倒问题

Notes on the Runge-Kutta method

示例：MMA求解经典力学问题

$U=-A/r^\alpha$ 下的经典粒子运动问题。保守力场下拉格朗日量为

\begin{matrix} (30) & L = {\dot{x}}^{2} (t) / 2 + {\dot{y}}^{2} / 2 - U \end{matrix}

拉格朗日方程

\begin{matrix} (31) & {\begin{cases} 0 = & \frac{d}{d t} (\frac{\partial L}{\partial \dot{x}}) - \frac{\partial L}{\partial x} \\ 0 = & \frac{d}{d t} (\frac{\partial L}{\partial \dot{y}}) - \frac{\partial L}{\partial y} \end{cases} \end{matrix}

$x(0),\dot{x}(0);y(0),\dot{y}(0)$ 。

1
1

Computational Physics

Episode 1. 计算方法

Computational Physics

Review：数值方法

插值

Lagrange插值

Newton-Cotes插值

其他插值方法

Runge-Kutta方法

求解二阶ODE

TEST

TEST

TEST