计算物理复习[ep.3|有限元方法、蒙特卡洛方法]

本文含有大量数学公式，请使用电脑或平板电脑浏览以获得最佳阅读体验

第六章 有限元方法

欧拉方程（欧拉-拉格朗日方程）

这是一种神奇的解法，可以把变分问题转化为微分方程来求解，让我们来看看这是如何实现的。

考虑一个最简形式泛函的变分问题：

$$\begin{aligned} \operatorname{min}I[y(x)]& =\int_{x_{1}}^{x_{2}}F{\big(}x,y,y'{\big)}\mathrm{d}x \\ M(I)& =\left\{y(x)\in C^1[x_0,x_1]\mid y(x_0)=y_0,y(x_1)=y_1\right\} \end{aligned} \tag{1}$$

若给函数 $y(x)$ 一个变分 $\delta y$ ，其泛函的变化量为：

$$\Delta I=I[y+\delta y]-I[y]=\int_{x_1}^{x_2}\big[F\big(x,y+\delta y,y'+\delta y'\big)-F\big(x,y,y'\big)\big]\mathrm{d}x \tag{2}$$

我们接下来将类比导数为零函数取极值，当泛函变化量在局域近似为0时，可认为到达局部最优值。

取原函数 $F{\big(}x,y,y'{\big)}$ 的梯度：

$$\nabla F\bigl(x,y,y'\bigr)=\left(\frac{\partial F}{\partial x},\frac{\partial F}{\partial y},\frac{\partial F}{\partial y'}\right)^\text{T} \tag{3}$$

由多元函数泰勒展开式 $f(\mathbf{x}_k+\mathbf{d}_k)\approx f(\mathbf{x}_k)+\mathbf{g}^T(\mathbf{x}_k)\mathbf{d}_k$，对 $F\bigl(x,y+\delta y,y^{\prime}+\delta y^{\prime}\bigr)$ 进行一阶泰勒近似：

$$F\bigl(x,y+\delta y,y'+\delta y'\bigr)=F\bigl(x,y,y'\bigr)+\left(\frac{\partial F}{\partial x},\frac{\partial F}{\partial y},\frac{\partial F}{\partial y'}\right)\cdot(0,\delta y,\delta y') \tag{4}$$

将（4）代入（2）得：

$$\Delta I=\int_{x_1}^{x_1}\bigg(\frac{\partial F}{\partial y}\delta y+\frac{\partial F}{\partial y'}\delta y'\bigg)\mathrm{d}x \tag{5}$$

考虑扰动函数的任意性，可令：

$$\delta y=\epsilon\eta(x),\quad\eta(x_0)=\eta(x_1)=0 \tag{6}$$

其中 $\eta(x)$ 是扰动函数，因为有 $y(x_0)=y_0,y(x_1)=y_1$ 的约束条件，$\eta(x)$ 在 $x_0$ 和 $x_1$ 处为零，即 $\delta y|_{x_0,x_1}=0$，$y(x)$ 在 $x_0$ 和 $x_1$ 处不能动。$\epsilon$ 是个放缩系数，在这里可以看作是扰动的幅度。（5）式因此变成：

$$\Delta I=\epsilon\int_{x_1}^{x_2}\bigg(\frac{\partial F}{\partial y}\eta+\frac{\partial F}{\partial y'}\eta'\bigg)\mathrm{d}x \tag{7}$$

你可不要小看这个 $\epsilon$ 哦，就是这个 $\epsilon$ 的引入让我们可以对（5）式使用费马引理，即令 $\Phi(\epsilon)=I[y+\epsilon\eta]$ 时，有 $\Phi'(\epsilon)|_{\epsilon=0}=0$，从而获得泛函的极值条件 $\Delta I=0$：

$$\int_{x_1}^{x_2}\bigg(\frac{\partial F}{\partial y}\eta+\frac{\partial F}{\partial y'}\eta'\bigg)\mathrm{d}x=0 \tag{8}$$

费马引理：设函数 $f(x)$ 在点 $x_0$ 的某邻域 $U(x_0)$ 内有定义，并且在 $x_0$ 处可导，如果对任意的 $x∈U(x_0)$，有 $f(x)\le f(x_0)\space or \space f(x)\ge f(x_0)$，那么 $f'(x_0)=0$

Cite From Wikipedia

对一阶导数 $\eta'$ 进行分部积分：

$$\int_{x_1}^{x_2}\bigg(\frac{\partial F}{\partial y'}\eta'\bigg)\mathrm{d}x=\frac{\partial F}{\partial y'}\eta\bigg|_{x_1}^{x_2}-\int_{x_2}^{x_2}\bigg(\frac{\mathrm{d}}{\mathrm{d}x}\frac{\partial F}{\partial y'}\eta\bigg)\mathrm{d}x \tag{9}$$

由于 $\eta(x_0)=\eta(x_1)=0$，等号右边第一项等于零，上式变为：

$$\int_{x_{1}}^{x_{2}}\eta\bigg(\frac{\partial F}{\partial y}-\frac{\mathrm{d}}{\mathrm{d}x}\frac{\partial F}{\partial y^{\prime}}\bigg)\mathrm{d}x=0 \tag{10}$$

因为 $\eta$ 是任意的，不一定为零，取括号内公式令其为零：

$$\frac{\partial F}{\partial y}-\frac{\mathrm{d}}{\mathrm{d}x}\frac{\partial F}{\partial y'}=0 \tag{11}$$

（11）式就是我们要求的欧拉（变分）方程啦！

微分方程边值问题转化为变分问题

既然一个变分问题可以被转化为微分方程问题，那么同样可以把一个微分问题转化为变分问题来求解，具体分以下两步走：

将微分方程转化为弱形式

以下述方程为例：

$$\begin{cases}-[p(x)y']'+q(x)y=f(x),x_0<x<x_1\\ y(x_0)=y_0,y(x_1)=y_1\end{cases} \tag{12}$$

两侧同时乘上测试函数（对变分问题来说是扰动函数 $\eta$ ） $z(x),s.t.\ z(x)|_{x_0,x_1}=0$

$$-[p(x)y']'z(x)+q(x)z(x)y=f(x)z(x) \tag{13}$$

等号两侧同时积分：

$$\int_{x_0}^{x_1}-\left[p(x)y'\right]'z(x)\mathrm dx+\int_{x_0}^{x_1}q(x)z(x)y\mathrm dx=\int_{x_0}^{x_1}f(x)z(x)\mathrm dx \tag{14}$$

对高阶导数项进行分部积分：

$$\int_{x_{0}}^{x_{1}}-\left[p(x)y^{\prime}\right]^{\prime}z(x)\mathrm{d}x=-p(x)y^{\prime}z(x)|_{x_{0}}^{x_{1}}+\int_{x_{0}}^{x_{1}}\left[p(x)y^{\prime}\right]z^{\prime}(x)\mathrm{d}x \tag{15}$$

回代到式（14）中并做移项：

$$\int_{x_0}^{x_1}\big\{p(x)y'(x)z'(x)+q(x)y(x)z(x)-f(x)z(x)\big\}\mathrm{d}x=0 \tag{16}$$

至此，原二阶微分方程初值问题就变为以下仅需一阶可导的弱形式问题：

在集合 $M(I)=\left\{y(x)\in C^1[x_0,x_1]\mid y(x_0)=y_0,y(x_1)=y_1\right\}$ 中找到一个函数 $y(x)$，使得对于任意测试函数 $z(x)$ 均有：

$$\int_{x_0}^{x_1}\big\{p(x)y'(x)z'(x)+q(x)y(x)z(x)-f(x)z(x)\big\}\mathrm{d}x=0$$

已经有一点变分问题的意思了对不对？

构造等价变分问题

由欧拉方程的推导可知：若令 $\Phi(\epsilon)=I[y+\epsilon\eta]$ ，当 $\Phi'(\epsilon)|_{\epsilon=0}=0$ 时泛函 $I[y(x)]$ 取得极小值，现构造泛函：

$$I[y(x)]=\int_{x_0}^{x_1}\big[p(x)y'^2(x)+q(x)y^2(x)-2f(x)y(x)\big]\mathrm{d}x \tag{17}$$

其满足：

$$\lim_{\epsilon \to 0} \Phi'(\epsilon) =\frac{I[y + \epsilon z] - I[y]}{\epsilon - 0} = \int^{x_1}_{x_0}[2p(x)y^{\prime}(x)z^{\prime}(x)+2q(x)y(x)z(x)-2f(x)z(x)]\mathrm{d}x \tag{18}$$

于是上述弱形式问题转化为以下变分问题：

$$\begin{aligned} \operatorname{min}I[y(x)]=\int_{x_{0}}^{x_{1}}\left[p(x)y^{\prime2}(x)+q(x)y^{2}(x)-2f(x)y(x)\right]\mathrm{d}x \\ M(I)=\left\{y(x)\in C^{1}[x_{0},x_{1}]\mid y(x_{0})=y_{0},y(x_{1})=y_{1}\right\} \end{aligned} \tag{18}$$

转化后的变分问题有以下特点：

• 变分问题只不过是原来偏微分方程做分部积分后的等价形式
• 只不过变分问题的解可以更弱。比如二阶方程，解应该二阶可导，对于变分问题，一阶可导就行。

Galerkin 方法

基函数系的选择

把微分问题转化为变分问题只是求解的第一步，如果想要得到目标函数 $y(x)$ ，我们还要接着解决变分问题。变分问题所需求解的函数为连续函数，不利于通过数值方法进行求解。在实际求解时，一般使用 Galerkin 方法，将用一组基函数的有限维线性组合近似表示待求连续函数。

以下面变分问题为例：

$$\begin{aligned} \operatorname{min}I[y(x)]=\int_{x_{0}}^{x_{1}}\left[p(x)y^{\prime2}(x)+q(x)y^{2}(x)-2f(x)y(x)\right]\mathrm{d}x \\ M(I)=\left\{y(x)\in C^{1}[x_{0},x_{1}]\mid y(x_{0})=y_{0},y(x_{1})=y_{1}\right\} \end{aligned}$$

要寻找一组基函数 $\varphi_k(x)$ 用其线性组合求解 $y(x)\approx\sum_{k=1}^na_k\varphi(x)$，这组基函数有以下性质：

• 每一个 $\varphi_k(x)$ 在 $[x_0,x_1]$ 上连续可导

• 每一个 $\varphi_k(x)$ 都满足边界条件 $\varphi_k(x_0)=\varphi_k(x_1)=0$

• 这组基函数 $\varphi_k(x)$ 中任意有限个函数互相线性无关

• 基函数 $\varphi_k(x)$ 的集合是完备的

• 对于任意给定的 $\varepsilon>0$ ，存在一个整数 $N$ 和常数组 $\{a_k\}^N_{k-1}$ ，使得：

  $$\begin{gathered} \operatorname*{max}_{x\in[x_{0},x_{1}]}\Big|\varphi(x)-\sum_{k=1}^{N}a_{k}\varphi_{k}(x)\Big|<\varepsilon \\ \operatorname*{max}_{x\in[x_{0},x_{1}]}\bigg|\varphi^{\prime}(x)-\sum_{k=1}^{N}a_{k}\varphi_{k}^{\prime}(x)\bigg|<\varepsilon \end{gathered}$$

本案例中选择多项式函数系：$\{\varphi_k(x)\}=\left\{(x_1-x)(x-x_0)^k\right\}_{k=1}^\infty$

看到这里，我不由得想起当时参加量子力学课程的时候，学到用一组基函数对波函数进行展开那一部分，从而引入狄拉克符号 $|\psi\rangle=\sum_ic_i|u_i\rangle$ 然后就是各种把我搞得死去活来的运算，但是不学量子力学的话我将很难理解 Galerkin 方法是如何work的。这种线性展开的方法用处非常广，如果把一个英文的句子看作一个函数，各个字母在单词乃至句子中的位置是展开系数的话，那（A-Z）英文字母表就是一组可以用来展开句子的互相线性无关（任意一个字母无法表示除它自己之外的其他字母）又完备的基函数，希望这能帮助你理解 Galerkin 方法到底干了些什么。

计算函数近似表达式

基函数已知，如果可以求出其对应的展开系数向量 $\mathbf{a}=(a_1,a_2,\cdots,a_n)^T$ ，就能得到线性组合的表达式。把线性组合代入泛函（17）：

$$I[y(x)]=\int_{x_q}^{x_1}\left[p(x)\left(\sum_{k=1}^{n}a_k\varphi'_k(x)\right)^2+q(x)\left(\sum_{k=1}^{n}a_k\varphi_k(x)\right)^2-2f\sum_{k=1}^{n}a_k\varphi_k(x)\right]\mathrm{d}x \tag{19}$$

显然，上式为一个关于 $a_k$ 的二次多项式，各项系数分别为：

$$\begin{aligned} a_{i}a_{j}&: \int_{x_0}^{x_1}\big[p(x)\varphi'_i(x)\varphi'_j(x)+q(x)\varphi_i(x)\varphi_j(x)\big]\mathrm{d}x \\ a_{i}&: -2\int_{x_0}^{x_1}f(x)\varphi_i(x)\mathrm{d}x \end{aligned} \tag{20}$$

令：

$$\begin{aligned} \alpha_{i,j} &=\int_{x_{0}}^{x_{1}}\big[p(x)\varphi_{i}^{\prime}(x)\varphi_{j}^{\prime}(x)+q(x)\varphi_{i}(x)\varphi_{j}(x)\big]\mathrm{d}x \\ \beta_{k} &=\int_{x_0}^{x_1}f(x)\varphi_k(x)\mathrm dx \end{aligned} \tag{21}$$

则（19）可表示为展开系数向量 $\mathbf{a}$ 的函数：

$$I[y_n(x)](\mathbf{a})=\sum_{i=1}^{n}\sum_{j=1}^{n}\alpha_{i,j}a_i a_j-2\sum_{k=1}^{n}\beta_k a_k \tag{22}$$

化为矩阵形式：

$$I[y_{n}(x)](\mathbf{a})=\mathbf{a}^{T}\mathbf{A}\mathbf{a}-2\mathbf{b}^{T}\mathbf{a} \tag{23}$$

其中：

$$\mathbf{A}={\left(\begin{array}{l l l l}{\alpha_{11}}&{\alpha_{12}}&{\cdots}&{\alpha_{1n}}\\ {\alpha_{21}}&{\alpha_{22}}&{\cdots}&{\alpha_{2n}}\\ {\vdots}&{\vdots}&{\ddots}&{\vdots}\\ {\alpha_{n1}}&{\alpha_{n2}}&{\cdots}&{\alpha_{n n}}\end{array}\right)}\quad\mathbf{b}={\left(\begin{array}{l}{\beta_{1}}\\ {\beta_{2}}\\ {\vdots}\\ {\beta_{n}}\end{array}\right)} \tag{24}$$

欲极小化该泛函，可直接求其梯度并令梯度为0：

$$\begin{array}{c} \nabla I(\mathbf a)=2\mathbf A\mathbf a-2\mathbf b=\mathbf0\to\mathbf A\mathbf a=\mathbf b \\ \Downarrow \\ \left(\begin{array}{l l l l}{\alpha_{11}}&{\alpha_{12}}&{\cdots}&{\alpha_{1n}}\\ {\alpha_{21}}&{\alpha_{22}}&{\cdots}&{\alpha_{2n}}\\ {\vdots}&{\vdots}&{\ddots}&{\vdots}\\ {\alpha_{n1}}&{\alpha_{n2}}&{\cdots}&{\alpha_{n n}}\end{array} \right)\left(\begin{array}{l}{a_{1}}\\ {a_{2}}\\ {\vdots}\\ {a_{n}}\end{array}\right)= \left(\begin{array}{l}{\beta_{1}}\\ {\beta_{2}}\\ {\vdots}\\ {\beta_{n}}\end{array}\right) \end{array}$$

至此，我们就通过把目标函数用一组基函数展开 $y(x)=\sum_{k=1}^na_k\varphi(x)$ 把一个复杂的变分问题转化为了简单的求解方程组的问题。这不能不说是很妙的，尤其的最后一步令泛函梯度为零，这一步直接将一个泛函转化为了矩阵之间的关系。Galerkin 方法为很多微分方程边值问题提供了一个简单的解决方案，当计算流程明晰之后，我们所要做的工作无非就是选择一组最佳的基函数罢了。

一维有限元方法

上一个小结里提到，所有用 Galerkin 方法实现的解微分方程边值问题的流程都十分相似，只是在基函数的选取上有些不同，有限元就是这样。我们将选择一组特殊的基函数对目标函数进行展开，它的表达式是这样的：

$$\varphi_i(x)=\begin{cases}\frac{(x-x_{i-1})}{h},&x_{i-1}< x\leqslant x_i\\ \frac{(x_{i+1}-x)}{h},&x_i< x < x_{i+1}\\ 0,&else \end{cases}(i=1,2,\cdots,N) \tag{25}$$

很明显这是个分段函数，它长这个模样：

不同的颜色对应着不同的 $i$ ，可以看到，这里我一共分了8个函数。下面我给出了一个有限元方法求解常微分方程边值问题的案例，例题相关代码在这里，该案例包括把微分方程转化为变分问题和有限元求解变分问题这一整个流程。

例题在这里！

使用一维有限元方法求解下列常微分方程边值问题

$$\left\{\begin{array}{l}y'' - y + x = 0, \ \ x\in[1,2]\\y(1) = 1,y(2) = 3\end{array}\right.\tag{1}$$

将微分方程转化为弱形式

对原函数 $y'' - y + x = 0$ 变形为 $-y'' + y = x$

两侧同时乘以测试函数 $z(x), s.t.\ z(1) = 0,z(2) = 0$

$$-y''z(x) + yz(x) = xz(x)\tag{2}$$

等式两侧同时积分：

$$\int^2_1-y''z(x)\mathrm{d}x + \int^2_1yz(x)\mathrm{d}x = \int^2_1xz(x)\mathrm{d}x\tag{3}$$

高阶导数项进行分部积分：

$$\int^2_1-y''z(x)\mathrm{d}x = -y'z(x)|^2_1 + \int^2_1y'z'(x)\mathrm{d}x\tag{4}$$

因为 $z(1) = 0,z(2) = 0$ ，上式中 $-y'z(x)|^2_1 = 0$

把（4）式代入（3），移项得：

$$\int^2_1\{y'z'(x) + yz(x) - xz(x)\}\mathrm{d}x = 0\tag{5}$$

则原微分方程初值问题变成以下弱形式问题：

在集合$M(I) = \{y(x)\in C^1[1,2]|y(1) = 1,y(2) = 3 \}$ 中找到一个函数 $y(x)$ ，使得对于任意测试函数 $z(x), s.t.\ z(1) = 0,z(2) = 0$ 均有：

$$\int^2_1\{y'z'(x) + yz(x) - xz(x)\}\mathrm{d}x = 0$$

构造等价变分问题

令$\Phi(\epsilon) = I[y + \epsilon z]$，有：

$$\Phi'(\epsilon) =\frac{I[y + \epsilon z] - I[y]}{\epsilon - 0}\tag{6}$$

其中，泛函 $I[y]$ 取极小值的条件为：

$$\Phi'(\epsilon)|_{\epsilon = 0} = 0\tag{7}$$

据此，我们构造泛函 $I[y(x)]$ :

$$I[y(x)] = \int^2_1[y'^2(x) + y^2(x) - 2xy(x)]\mathrm{d}x\tag{8}$$

其满足：

$$\lim_{\epsilon \to 0} \Phi'(\epsilon) =\frac{I[y + \epsilon z] - I[y]}{\epsilon - 0} = \int^2_1[2y'z' + 2yz - 2xz]\mathrm{d}x\tag{9}$$

于是上述弱形式问题转化为求解以下变分问题：

$$\begin{array}{c}\mathrm{min}\ I[y(x)] = \int^2_1[y'^2(x) + y^2(x) - 2xy(x)]\mathrm{d}x\\M(I) = \{y(x)\in C^1[1,2]|y(1) = 1,y(2) = 3 \}\end{array}\tag{10}$$

有限元求解

对于要求解的线性函数：

$$\mathbf{Aa} = \mathbf{b}\tag{11}$$

由基函数的形状可得，当 $|i - j| > 1$ 时，矩阵元 $\alpha_{i,j}$ :

$$\alpha_{i,j} = \int^2_1[\varphi'_i(x)\varphi'_j(x) + \varphi_i(x)\varphi_j(x)]\mathrm{d}x = 0 \tag{12}$$

则系数矩阵 $\mathbf{A}$ 为三对角矩阵：

$$\mathbf{A} = \begin{bmatrix} \alpha_{1,1} & \alpha_{1,2} & 0 &\cdots & 0 \\ \alpha_{2,1} & \alpha_{2,2} & \alpha_{2,3} &\cdots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & \alpha_{n-1,n-2} & \alpha_{n-1,n-1} & \alpha_{n-1,n} \\ 0 & \cdots & 0 & \alpha_{n,n-1} & \alpha_{n,n} \end{bmatrix}\tag{13}$$

向量 $\mathbf{b}$ 为：

$$\mathbf{b} = \begin{bmatrix}\beta_{1} \\\beta_{2} \\\vdots \\\beta_{n - 1} \\\beta_{n}\end{bmatrix}\tag{14}$$

由于只有 $\varphi_1(x)$ 对边界条件 $y(1) = 1$ 有贡献，则 $\mathbf{A}$ 的第一行和 $\mathbf{a}$ 相乘后应该为 $\mathbf{b}$ 向量的第一个分量，该分量应为1。为满足此条件，矩阵 $\mathbf{A}$ 的第一行应该变为除了 $\alpha_{1,1}=1$ 之外皆为 0 的状态。此处令 $\alpha_{1,1}=1$ 是为了求解 $\mathbf{a}$ 向量真正的数值。

同理，考虑另一个边界条件 $y(2) = 3$ ，矩阵 $\mathbf{A}$ 的第 n 行应该变为除了 $\alpha_{n,n}=1$ 之外皆为 0 的状态，即：

$$\mathbf{A} = \begin{bmatrix} 1 & 0 & 0 &\cdots & 0 \\ \alpha_{2,1} & \alpha_{2,2} & \alpha_{2,3} &\cdots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & \alpha_{n-1,n-2} & \alpha_{n-1,n-1} & \alpha_{n-1,n} \\ 0 & \cdots & 0 & 0 & 1\end{bmatrix}\tag{14}$$

此时向量 $\mathbf{b}$ 为：

$$\mathbf{b} = \begin{bmatrix}1 \\\beta_{2} \\\vdots \\\beta_{n - 1} \\3\end{bmatrix}\tag{15}$$

剩下步骤利用追赶法快速求解即可。

第七章 蒙特卡洛方法

随机数采样方法：逆变换法与舍取法

已知随机变量 $X$ 的概率密度函数为 $f(x)$，利用 $[0,1]$ 均匀随机数生成分布满足 $f(x)$ 的随机数。

$$\xi\sim\mathrm{U}[0,1]\to\eta\sim f(x) \tag{26}$$

以下是两种用于取样的方法：

逆变换法

利用 $f(x)$ 求出 $\eta$ 的累积分布函数，可以证明：

$$F(\eta)=\int_{-\infty}^\eta f(x)\mathrm dx\sim\mathrm U[0,1] \tag{27}$$

$$F(\eta)=\int_{-\infty}^\eta f(x)\mathrm{d}x=\xi\quad\Longrightarrow\quad\eta=F^{-1}(\xi) \tag{28}$$

下面我们用一个例子来说明如何使用逆变换法：

利用 $[0,1]$ 均匀随机数生成指数分布的随机变量

$$f(x)=\begin{cases}\lambda^{-1}\mathrm e^{-x/\lambda},&x>0,\lambda>0\\ 0,&x\leqslant0\end{cases}$$

先根据概率密度函数求出累积分布函数：

$$F(\eta)=\int_{-\infty}^{\eta}f(t)\mathrm{d}t=\int_{0}^{\eta}\lambda^{-1}\mathrm{e}^{-t/\lambda}\mathrm{d}t=1-\mathrm{e}^{-\eta/\lambda}=\xi$$

求反函数：

$$F^{-1}(\xi) =-\lambda\ln(1-\xi)=\eta$$

示例代码和结果：

% 用[0,1]均匀分布获得指数函数
n = 10000;
lambda = 0.5;
X = zeros(1,n);
for i = 1:n
    X(i) = -1*lambda*log(1-rand()); % 反函数
end
% 以下内容↓只和画图有关
h = histogram(X);
txt = ['$f(x)=' num2str(lambda) '^{-1}e^{-x/' num2str(lambda) '}$'];
xlim([0 max(X)])
text(0.6*max(X),0.8*max(h.Values),txt,'Interpreter','latex','FontSize',14);
title("Exponential Distribution")

舍选法

当无法给出累积分布函数逆函数的解析表达式时，舍选法是另外的选择。舍选法的适用范围比逆变换法要大，只要给出概率密度函数（PDF）的解析表达式即可，而大多数常用分布的PDF是可以查到的。

算法流程：

1. 确定概率密度函数 $f(x)$ 的区间 $[x_{min},x_{max}]$ 及其最大最小值 $[y_{min},y_{max}]$
2. 生成一个均匀分布随机数：$x\sim U[x_{min},x_{max}]$
3. 独立的生成另一个均匀分布随机数：$y\sim U[y_{min},y_{max}]$
4. 如果 $y \le f(x)$ ，保留生成的点，否测重复上述步骤。

如图，绿色的点在 PDF 下面，接受；红色的点在 PDF 下面，舍去。

示例代码和结果：

%舍选法求任意函数的分布(x)
clear, clc
n = 100000;
h = 0.01;
x = 0:h:1;
f = @(x)(20*x.*(1 - x).^3); % 概率密度函数（PDF）
y = f(x);
ymax = max(y);
ymin = min(y);

vec = [];
for i = 1:n
    yrand = ymin + (ymax - ymin)*rand; % 逆变换法求均匀分布
    R = rand; 
    if yrand <= f(R)
        vec = [vec R]; % 小于f(x)就保存
    end
end
% 以下内容↓只和画图有关
hold on
yyaxis left
histogram(vec,'FaceAlpha',0.2,'FaceColor','cyan')
yyaxis right
plot(x,y,'r-','LineWidth',1.5)
box on
title("Acceptance-Rejection Method")
text(0.6,0.8*ymax,"$f(x)=20x(1-x)^3$","Interpreter","latex",'FontSize',14)
legend("Distribution", "PDF")

蒙特卡洛求解方程、定积分与微分方程

求解方程

蒙特卡洛求根的思想是把“求根”这个问题转化为求概率的问题，设求：

$$f(x)=0\quad x\in[a,b] \tag{29}$$

设 $x^*$ 为该方程在 $[a,b]$ 上的解，则 $x^*$ 左边的区间 $[a,x^*]$ 长度为 $x^*-a$，若在 $[a,b]$ 上均匀取点（$\xi\sim U[a,b]$），则该点落到左边区间的概率为：

$$p=\frac{x^*-a}{b-a} \tag{30}$$

以此建立方程的根 $x^*$ 和概率的关系：

$$x^*=a+(b-a)p \tag{31}$$

若在 $[a, b]$ 区间上均匀取点 $N$ 次，有 $M$ 次满足 $f(x)\cdot f(a)>0$ （同一侧）：

$$\begin{array}{c} p\approx\frac{M}{N} \\ \Downarrow \\ x^*\approx a+(b-a)\frac{M}{N} \end{array} \tag{32}$$

求定积分（投点法 – 适用于 f (x) 恒正/负）

计算流程：

1. 找到 $f(x)$ 最大值 $f(x_m)$

2. 生成两个独立的 $U[0,1]$ 随机变量 $\xi_i,\eta_i$

3. 投点：

   $$\begin{cases}x_i=a+(b-a)\xi_i\\ y_i=f(x_m)\eta_i\end{cases} \tag{33}$$

4. 如果 $y_i<f(x_i)$：$M=M+1$，积分结果：

   $$I\approx\frac{M}{N}S=\frac{M}{N}(b-a)f(x_m) \tag{34}$$

这个方法看起来很复杂，其实很好理解。想象你有一个矩形的盘子（上图）面积为 $S$（这里的 $S=1\times0.7$），盘子的左下角在积分的下限，右下角在积分的上限，你在盘子里画了一条红色的线，是你将要积分的函数的函数图像。你开始向盘子里扔豆子，在我们的积分面积（本案例中是红线下方的区域）里的豆子标记为绿色，外面的是蓝色。假设你一共扔了 $N$ 颗豆子，其中有 $M$ 颗是绿色，那么豆子标记为绿色的概率和盘子面积的乘积就可以近似为积分的结果，即 $I\approx\frac{M}{N}S$。这种方法虽然有诸多限制，但是这个方法也有一些优点，特别是用于计算一些形状不规则的 封闭 图形的面积。只要数清楚落在这个图形里豆子的数目，除以扔的所有豆子的数目再乘上总撒豆面积就是这个不规则图形面积的近似结果了。

求解微分方程

用蒙特卡洛求解微分方程时，实质上是一个复化求积的过程，只不过在子区间是用算术平均方法来求积。关于复化求积可以看我的这篇文章。

设：

$$\frac{\mathrm{d}y}{\mathrm{d}x}=f(x,y)\quad y(x_0)=y_0 \tag{35}$$

求积分：

$$y(x) = y(x_0) + \int_{x_0}^xf(t,y)\mathrm{d}t \tag{36}$$

对待求区间进行剖分：

$$y(x)=y(x_0)+\sum\limits_{i=1}^N\int_{x_{i-1}}^{x_i}f(t,y)\mathrm{d}t \tag{37}$$

接下来对子区间进行积分，设 $u(x) = \frac{1}{x_i-x_{i-1}}$ 为 $\xi\sim U[x_i,x_{i-1}]$ 的概率密度函数，则：

$$I_i = \int_{x_{i-1}}^{x_i}f(t,y)\mathrm{d}t = \int_{x_{i-1}}^{x_i}\frac{f(t,y)}{u(t)}u(t)\mathrm{d}t \tag{38}$$

令：

$$g(t,y) =\frac{f(t,y)}{u(t)} = (x_{i}-x_{i-1})f(t,y) \tag{39}$$

那么求子区间积分 $I_i$ 就变成了求 $g(t,y)$ 的期望：

$$I_i =\int_{x_{i-1}}^{x_i}g(t,y)u(t)\mathrm{d}t\approx\frac{1}{K}\sum_{k=1}^{K}g(t_k,y)=\frac{x_{i}-x_{i-1}}{K}\sum_{k=1}^{K}f(t_k,y) \tag{40}$$

最后，化为递推式：

$$y(x_i)=y(x_{i-1})+\frac{x_{i}-x_{i-1}}{K}\sum_{k=1}^{K}f(t_k,y(x_{i-1})) \tag{41}$$

以下是一个关于蒙特卡洛求解微分方程的例题，希望能让你更好地理解这个方法：

例题在这里！

蒙特卡洛方法求解以下一阶微分方程问题：

$$\left\{\begin{array}{l}\frac{\mathrm{d}y}{\mathrm{d}x}=x\sqrt{|y|}+(sin(\frac{\pi x}{2}))^3-3\Theta(x-2)\quad(-4<x<5)\\y(-4) = 4\end{array}\right .$$

其中 $\Theta(x)$ 为阶跃函数：

$$\Theta(x) = \begin{cases}1\quad x>0 \\0\quad else\end{cases}$$

计算流程和（35-41）相同，示例代码如下：

%蒙特卡洛方法解微分方程
clear, clc
a = -4;
b = 5;
h = 0.02;
x = a:h:b;
N = size(x, 2);
y = zeros(1, N);
y(1) = 4;
K = 10000;

for i = 2:N
    sum = 0;
    for k = 1:K
        xk = x(i - 1) + (x(i) - x(i - 1))*rand(); % x_{i-1}到x_i的均匀分布
        sum = sum + f(xk, y(i - 1));
    end
    sum = sum/K*(x(i) - x(i - 1));
    y(i) = y(i - 1) + sum; % 递推式
end

plot(x,y)

function z = f(x, y)
    z = x*sqrt(abs(y)) + sin(x*pi/2)^3 - 3*Theta(x - 2);
end

function Z = Theta(x) % 阶跃函数
    if x > 0
        Z = 1;
    else
        Z = 0;
    end
end

終了

下一篇是最后一篇了，因为后三章主要考概念所以公式不会像这几章这么多了。。。吧？坚持一下，就要结束了 TAT

有误之处欢迎反馈！