CS205A Lecture 18 and 19 PDE

课程主页：https://graphics.stanford.edu/courses/cs205a-13-fall/schedule.html

这次回顾最后两讲，这部分简单介绍了PDE。

Chapter 14 偏微分方程

这一讲介绍PDE，即偏微分方程，为了方便后续讨论，这里给出如下记号，假设

$f : \mathbb{R}^{3} \rightarrow \mathbb{R}, \vec{v} : \mathbb{R}^{3} \rightarrow \mathbb{R}^{3}$

那么

$\begin{aligned} \text{Gradient: } &\nabla f \equiv\left(\frac{\partial f}{\partial x_{1}}, \frac{\partial f}{\partial x_{2}}, \frac{\partial f}{\partial x_{3}}\right) \\ \text { Divergence: } & \nabla \cdot \vec{v} \equiv \frac{\partial v_{1}}{\partial x_{1}}+\frac{\partial v_{2}}{\partial x_{2}}+\frac{\partial v_{3}}{\partial x_{3}} \\ \text { Curl: }& \nabla \times \vec{v} \equiv\left(\frac{\partial v_{3}}{\partial x_{2}}-\frac{\partial v_{2}}{\partial x_{3}}, \frac{\partial v_{1}}{\partial x_{3}}-\frac{\partial v_{3}}{\partial x_{1}}, \frac{\partial v_{2}}{\partial x_{1}}-\frac{\partial v_{1}}{\partial x_{2}}\right) \\ \text { Laplacian: } & \nabla^{2} f \equiv \frac{\partial^{2} f}{\partial x_{1}^{2}}+\frac{\partial^{2} f}{\partial x_{2}^{2}}+\frac{\partial^{2} f}{\partial x_{3}^{2}} \end{aligned}$

高维情形的记号类似。

动机

继续讨论之前来看一个例子：

例1

定义

$E[f]=\int_{\Omega}\|\nabla f(\vec{x})\|_{2}^{2} d \vec{x}$

我们的目的是

$\begin{aligned} \text { minimize: }& E[f] \\ \text { such that: }& f(\vec{x})=g(\vec{x}) &\forall x \in \partial \Omega \end{aligned}$

如果$f$最小化$E$，那么对任意$h$，我们有

$E[f+h]\ge E[f]$

现在考虑$E[f+\epsilon h]$

$\begin{aligned} E[f+\epsilon h] &=\int_{\Omega}\|\nabla f(\vec{x})+\epsilon \nabla h(\vec{x})\|_{2}^{2} d \vec{x} \\ &=\int_{\Omega}\left(\|\nabla f(\vec{x})\|_{2}^{2}+2 \epsilon \nabla f(\vec{x}) \cdot \nabla h(\vec{x})+\epsilon^{2}\|\nabla h(\vec{x})\|_{2}^{2}\right) d \vec{x} \end{aligned}$

关于$\epsilon $微分可得

$\begin{aligned} \frac{d}{d \epsilon} E[f+\epsilon h] &=\int_{\Omega}\left(2 \nabla f(\vec{x}) \cdot \nabla h(\vec{x})+2 \epsilon\|\nabla h(\vec{x})\|_{2}^{2}\right) d \vec{x} \\ \frac{d}{d \epsilon} E\left.[f+\epsilon h]\right|_{\epsilon=0} &=2 \int_{\Omega}[\nabla f(\vec{x}) \cdot \nabla h(\vec{x})] d \vec{x} \end{aligned}$

上述推导对任意$h$都成立，特别的，我们取$h(\vec x)=0,\vec x \in \partial \Omega$，那么利用分步积分可得

$\frac{d}{d \epsilon} E\left.[f+\epsilon h]\right|_{\epsilon=0}=-2 \int_{\Omega} h(\vec{x}) \nabla^{2} f(\vec{x}) d \vec{x}$

上式应该恒等于$0$，所以

$\nabla^2 f(\vec x)= 0 , x\in \Omega\setminus \partial \Omega$

所以我们需要求解的PDE为

$\begin{aligned} \nabla^{2} f(\vec{x}) &=0 \\ f(\vec{x}) &=g(\vec{x}) &\forall \vec{x} \in \partial \Omega \end{aligned}$

该方程称为拉普拉斯方程。

基本记号

这里假设

$f : \mathbb{R}^{n} \rightarrow \mathbb{R}$

定义如下记号（以此类推）

$\begin{aligned} f_{x} & \equiv \frac{\partial f}{\partial x} \\ f_{y} & \equiv \frac{\partial f}{\partial y} \\ f_{x y} & \equiv \frac{\partial^{2} f}{\partial x \partial y} \end{aligned}$

根据边界条件分为以下三类：

Dirichlet条件给出$\partial \Omega$上的$f(\vec x)$
Neumann条件给出$\partial \Omega$上的$\nabla f(\vec x)$
Mixed或Robin条件给出上述两者的结合

对方程建模

实际中讨论最多的是二阶PDE

$\sum_{i j} a_{i j} \frac{\partial f}{\partial x_{i} \partial x_{j}}+\sum_{i} b_{i} \frac{\partial f}{\partial x_{i}}+c=0$

定义梯度算子

$\nabla \equiv\left(\frac{\partial}{\partial x_{1}}, \frac{\partial}{\partial x_{2}}, \ldots, \frac{\partial}{\partial x_{n}}\right)$

所以上述方程可以改写为

$\left(\nabla^{\top} A \nabla+\nabla \cdot \vec{b}+c\right) f=0$

如果$A$是正定或者负定矩阵，那么该系统是椭圆型。
如果$A$是半正定或者半负定矩阵，那么该系统是抛物型。
如果$A$只存在一个特征值和其他特征值符号不同，那么该系统是双曲型。
如果不满足上述条件，那么称该系统是超双曲型。

椭圆型PDE

考虑拉普拉斯方程：

$\nabla^2 f=g$

抛物型PDE

$\frac{\partial f}{\partial t}=\alpha \nabla^{2} f$

双曲型PDE

$\frac{\partial^{2} f}{\partial t^{2}}-c^{2} \nabla^{2} f=0$

将微分看成算子

由上一讲的思路，可以得到$f(x)$的二阶导数的近似公式

$f^{\prime \prime}(x)=\frac{f(x+h)-2 f(x)+f(x-h)}{h^{2}}+O(h)$

现在假设$f(x)$在$[0,1]$上有样本

$y_{0} = f(0), y_{1}= f(h), y_{2}= f(2 h), \ldots, y_{n}=f(n h)$

那么利用上述公式可以得到

$y_{k}^{\prime \prime} \approx \frac{y_{k+1}-2 y_{k}+y_{k-1}}{h^{2}}$

另一方面，考虑二元函数$f :[0,1] \times[0,1] \rightarrow \mathbb{R}$，可以用类似的方法得到

$\left(\nabla^{2} y\right)_{k, \ell} =\frac{1}{h^{2}}\left(y_{(k-1), \ell}+y_{k,(\ell-1)}+y_{(k+1), \ell}+y_{k,(\ell+1)}-4 y_{k, \ell}\right)$

求解

将之前的微分算子转化为矩阵，然后求解线性方程组即可，这里从略。