最近在做的项目是关于偏微分方程数值解的快速计算，于是先学一下偏微分方程的基本原理。这篇文章主要是对于3b1b的微分方程系列的简要总结

常微分方程

案例

对问题的讨论起始于经典的单摆问题

image-20250111155903222

单摆角度θ\theta以右为正，根据上面的力的分解可以得到

a=x¨=−gsin⁡θa = \ddot x=-g\sin\theta

而

x=θLx = \theta L

于是

θ¨=−gLsin⁡θ\ddot\theta=-\frac gL\sin\theta

加入阻力项，得到

θ¨=−μθ˙−gLsin⁡θ\ddot\theta=-\mu\dot\theta-\frac gL\sin\theta

这就得到了一个描述单摆运动的二阶常微分方程（ODE）

这里我们用上方加点x˙\dot x来表示对时间求导ddtx\frac {\rm d}{\mathrm{d} t}x，也就是，使用物理中常用的记号

对案例的分析

视频中使用了相空间（Phase Space）来讨论这个问题，具体地讲，相空间是当前系统状态的空间，也就是，相空间中的每个点对应唯一确定的系统状态

在上面的单摆问题中，系统状态可以由θ\theta和θ˙\dot\theta唯一确定，于是相空间就是(θ,θ˙)(\theta,\dot\theta)组成的二维空间

为了使用相空间分析ODE，我们需要计算相空间中状态随时间的变化，也就是

ddt[θθ˙]=[θ˙θ¨]=[θ˙−μθ˙−gLsin⁡θ]\frac{\rm d}{\mathrm{d}t}\begin{bmatrix}\theta\\\dot\theta\end{bmatrix}= \begin{bmatrix}\dot\theta\\\ddot\theta\end{bmatrix}=\begin{bmatrix}\dot\theta\\-\mu\dot\theta-\frac gL\sin\theta\end{bmatrix}

于是，我们将系统状态随时间的变化表达为了系统状态自身的函数，也就是，状态的导数是一个向量场，在每一个确定的状态点，系统有确定的走向

事实上，如果把θ\theta和θ˙\dot\theta看成两个独立的变量，上面的做法相当于把二阶常微分方程转化为了二元一阶常微分方程组

image-20250111160912088

进行上面向量场分析的意义是，我们可以将向量场以图形的方式描述出来，从而通过给定起始点在相空间中的运动路径，来分析系统的变化

热传导方程

接下来让我们进入**偏微分方程（PDE）**的领域，讨论最简单也最常被讨论的偏微分方程——热传导方程

考虑一根初始温度不均匀但材料均匀的一维杆，随着时间推移，其温度会逐渐趋于均匀，于是杆上的温度是杆上的位置和时间两个变量的函数：

T(x,t)T(x,t)

接下来我们考虑热传导中每个固定点的温度变化，也就是，我们希望了解∂∂tT(x,t)\frac{\partial}{\partial t}T(x,t)与什么有关？

首先考虑：满足什么条件的点的温度会保持不变，即上面偏导数为零？

我们首先考虑离散的情形，假设杆上有NN个点，其中每个点受到的来自两侧两个点的热量应该是正比于它与两侧的点的温度差的，即：

ΔTkleft=−a(Tk−Tk−1)\Delta T_k^{\rm left}=-a(T_k-T_{k-1})

ΔTkright=−a(Tk−Tk+1)\Delta T_k^{\rm right}=-a(T_k-T_{k+1})

于是

ΔTk=−a(2Tk−Tk−1−Tk+1)=a[(Tk+1−Tk)−(Tk−Tk−1)]\Delta T_k=-a(2T_k-T_{k-1}-T_{k+1})=a[(T_{k+1}-T_k)-(T_k-T_{k-1})]

这实际上是TT在空间上的二阶差分，即差分的差分

类推到连续的情形，二阶差分将变成二阶微分，即温度对时间的偏导数与温度对位置的二阶偏导数成正比，这就得到了一维热传导方程：

∂T∂t(x,t)=α∂2T∂x2(x,t)\frac{\partial T}{\partial t}(x,t)=\alpha \frac{\partial^2 T}{\partial x^2}(x,t)

这里的二阶导数实际表示了一个点与它的近邻点的均值的差距，或者这个点相比它的近邻点的突出程度，这样理解就符合热传导的直观了

多维情形的推广非常简单：

∂T∂t=α(∂2T∂x2+∂2T∂y2+∂2T∂z2)=α∇2T\frac{\partial T}{\partial t}=\alpha\Big( \frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}+\frac{\partial^2 T}{\partial z^2}\Big)=\alpha\nabla^2T

其中Laplace算子定义为

∇2=∂2∂x2+∂2∂y2+∂2∂z2\nabla^2=\frac{\partial^2 }{\partial x^2}+\frac{\partial^2 }{\partial y^2}+\frac{\partial^2 }{\partial z^2}

注意，它实际上是以下Nabla算子的衍生：

∇=[∂∂x∂∂y∂∂z]\nabla=\begin{bmatrix}\frac\partial{\partial x}\\\frac\partial{\partial y}\\\frac\partial{\partial z}\end{bmatrix}

∇2=∇T∇=∥∇∥2\nabla^2=\nabla^T\nabla=\|\nabla\|^2

Fourier对热传导方程的求解

对热传导方程的求解正是Fourier于19世纪上半叶提出的，这也是Fourier级数的重要应用之一

首先，考虑以下初始温度分布：

T(x,0)=sin⁡xT(x,0)=\sin x

显然

∂T∂t(x,0)=α∂2T∂x2(x,0)=−αsin⁡x\frac{\partial T}{\partial t}(x,0)=\alpha \frac{\partial^2 T}{\partial x^2}(x,0)=-\alpha\sin x

于是，经过微小时间间隔后，每一点的温度经过了等比例的缩小，于是其缩小的结果仍然是正弦函数，于是下一时刻的温度仍然是等比例的缩小

因此，可以推知

T(x,t)=f(t)sin⁡xT(x,t)=f(t)\sin x

代入原方程，有

f′(t)sin⁡x=−αf(t)sin⁡xf'(t)\sin x=-\alpha f(t)\sin x

即

f(t)=e−αtf(t)=e^{-\alpha t}

于是

T(x,t)=e−αtsin⁡xT(x,t)=e^{-\alpha t}\sin x

是一个合理解

显然，如果我们推广初始温度分布为任意的Fourier三角基函数，都可以解出形式很好的解析解，而这里的PDE是线性的，即：若两个函数T1T_1和T2T_2是热传导方程的解，那么其任意线性组合

λ1T1+λ2T2\lambda_1T_1+\lambda_2T_2

也是热传导方程的解，即PDE的解空间是线性空间

PDE不总是线性的。这里的PDE是线性的是由于方程左侧和右侧的算子都是线性的，很显然，只要在右侧的某一个TT外面套上sin⁡\sin，它就不再线性

于是，我们可以对初始温度分布求解Fourier级数，从而利用PDE解空间的线性，求得其解

但，事情没有这么简单

边界条件

事实上，如果只考虑PDE本身，那么

T(x,t)=cx, c≠0T(x,t)=cx,\ c\neq0

也是合理的解，因为其对xx的二阶偏导永远是00

但显然，这样的温度不均的系统是无法稳定的，究其原因是：虽然杆中间的点在温度等于左右两个点的均值时就是平衡的，但对于杆两端的点，它们有一侧没有点，因此需要重新考虑此时的平衡条件，即边界条件

如果认为边界点与外界没有热量交换，那么边界条件为

∂T∂t(0,t)=∂T∂t(L,t)=0, ∀t>0\frac{\partial T}{\partial t}(0,t)=\frac{\partial T}{\partial t}(L,t)=0,\ \forall t>0

为了使Fourier基函数满足上面的边界条件，实际上只需要选取

cos⁡nπxL,n=0,1,2,⋯\cos\frac{n\pi x}{L},\quad n=0,1,2,\cdots

作为Fourier基函数即可，得到的解簇为

e−αn2π2tL2cos⁡nπxLe^\frac{-\alpha n^2\pi^2 t}{L^2}\cos\frac{n\pi x}{L}

将上面的解簇线性组合，使得初始温度分布是给出的温度分布的Fourier级数，即完成了求解

总结

对于偏微分方程的问题，PDE本身往往是不够的，决定解的条件还有边界条件和初值条件