关于PDE,我的问题是缺乏一种系统的把握。我所有的训练都是在数理方程的课上得到的,这些训练只告诉我某些十分特定的方程的解是某些特殊函数,至于这些方程的特定究竟特定到了什么程度,具体是怎样的条件,虽然也讲,但并不重视,主要重视的是造成这些方程的物理问题。所以就形成一个对物理问题的解的感觉。但是这个感觉现在也已经消退了。所以要重新从这些基础的东西学起。
分类
物理上主要是考虑二阶PDE,比如说双曲型的波动方程,椭圆形的泊松方程,抛物线型的导热方程。然后边界条件也主要见到三类,dirichlet给定边界的函数值,neumann条件给定边界的导数,cauchy两者都给定。实际上解决问题的时候,有唯一解的必要条件对于椭圆和抛物型方程来说是边界条件有d和n的其中之一,而对双曲型,则得是c。
Cauchy data
一维最简单的情况,知道位置和速度,就能推演以后的运动。这个时候这两个信息就是cauchy data。二维的时候
我们给一个曲面,上面的函数值和导数是知道的,原则上
∂n∂ti∂2φ∂ti∂tj∂2φ= def nμtiν∂xμ∂xν∂2φ= def tiνtjν∂xμ∂xν∂2φ
是知道的,但就这些还求不出jacobian的全部信息,还得要知道
∂n∂n∂2φ=defnμnν∂xμ∂xν∂2φ
才行。求这个东西的trick是取∂xμ∂xν∂2φ=ϕ0μν+nμnνΦ,其中ϕ是满足前两组方程的解,Φ是这个解之外的“正交”的部分。对于方程
aμν(x)∂xμ∂xν∂2φ+( lower orders )=0
我们就知道
aμν(xi)∂xμ∂xν∂2φ+( known lower orders )=0
aμνnμnνΦ+(knownstuff)=0
如果aμνnμnν是0的话就坏了,我们就求不出来Φ了。对于这种的点我们叫characteristic (surface),因为它老是连成面。以后我就管他叫char。char对解方程其实是有好处的,不是完全的坏事。关于解的信息可以在char上面传递。看一个例子:
a(x,y)∂x∂u+b(x,y)∂y∂u+c(x,y)u=f(x,y)
它是(v⋅∇)u+cu=f式的方程,给出一个u的流。每条流线上dtdu+c(t)u(t)=f(t),其中t是使得dtdx=a(x,y),dtdy=b(x,y)的的参数化。给一条线Γ上的cauchy data,这个线和流线(也就是char)相交的点处的信息全知道,所以每条char上面的结果都能解出来。如果相切了,那多半就是无解。
波动方程
一、
波动方程是一个二阶的双曲方程。这种方程高次项可以因式分解
a∂x2∂2+2b∂x∂y∂2+c∂y2∂2=(α∂x∂+β∂y∂)(γ∂x∂+δ∂y∂)+ lower, =(γ∂x∂+δ∂y∂)(α∂x∂+β∂y∂)+ lower.
然后写成
(α∂x∂+β∂y∂)U1+F1=0
(γ∂x∂+δ∂y∂)U2+F2=0
当做一维的PDE来求解。波动方程就可以这样解。还有一个简单的办法是d’Alambert弄出来的,令
ξ=x+ct
η=x−ct
就知道
∂ξ∂=∂ξ∂x∂x∂+∂ξ∂t∂t∂=21(∂x∂+c1∂t∂)
∂η∂=21(∂x∂−c1∂t∂)
(∂x2∂2−c21∂t2∂2)=(∂x∂+c1∂t∂)(∂x∂−c1∂t∂)=4∂ξ∂η∂2=0
所以说
φ(x,t)=f(ξ)+g(η)=f(x+ct)+g(x−ct)。
现在我们有初值φ(x,0)≡φ0(x) and φ˙(x,0)≡v0(x),这个就是cauchy data。我们就知道
f(x)+g(x)c(f′(x)−g′(x))=φ0(x)=v0(x)
以下这个积分的过程相当于是把解的信息在char上propagate:
f(x)−g(x)=c1∫0xv0(ξ)dξ+A
然后就一路得到结果。
二、
还有一种办法可以展示积分变换怎样解这个问题:做Fourier变换
φ(x,t)=∫−∞∞2πdk{a(k)eikx−iωkt+a∗(k)e−ikx+iωkt}
边界条件
Φ(k)= def ∫−∞∞φ(x,t=0)e−ikxdx
χ(k)= def ∫−∞∞φ˙(x,t=0)e−ikxdx
对应
\Phi(k)=a(k)+a^\*(-k)
\chi(k)=i \omega_{k}\left(a^{\*}(-k)-a(k)\right)
然后就知道结果了。
三、
格林函数作为另外一种具备普遍性的办法,可以解这个问题。对于有源的方程
c21∂t2∂2φ−∂x2∂2φ=q(x,t)
我们解
(c21∂t2∂2−∂x2∂2)G(x,t;ξ,τ)=δ(x−ξ)δ(t−τ)
考虑到因果性,我们是在t上面积分,因为t<τ这一段已经知道是0了。
G(x,τ+ε;ξ,τ)dtdG(x,τ+ε;ξ,τ)=0=c2δ(x−ξ)
这就是τ右边这一段的cauchy data。从f(x+ct)+g(x−ct)形式的解代入
f(x)+g(x)c(f′(x)−g′(x))=φ0(x)=v0(x)
得出来就是
G(x,t;ξ,τ)=θ(t−τ)2c∫x−c(t−τ)x+c(t−τ)δ(ζ−ξ)dζ=2cθ(t−τ){θ(x−ξ+c(t−τ))−θ(x−ξ−c(t−τ))}
一个算例:不同维度上的爆炸
这是一个挺有意思的事情。考虑三维的一个声波
∂x2∂2ϕ+∂y2∂2ϕ+∂z2∂2ϕ−c21∂t2∂2ϕ=0
把方程写成流方程的形式,就知道一些物理量
v1ρ1P1=∇ϕ=−c2ρ0ϕ˙=c2ρ1
只考虑球面波。直接写
ϕ(r,t)=r1f(t−cr)+r1g(t+cr)
后一项扔掉,只考虑膨胀。记q˙是体积膨胀的速度,v(r,t)=4πr2q˙(t)。那么就有
4πr2q˙(t)=v1(r,t)=∂r∂ϕ=−r21f(t−cr)−rc1f′(t−cr)
考虑以下两种情况:
- 近场有−4π1q˙(t)=f(t)
- 远场有v1=∂r∂ϕ≈−rc1f′(t−cr)
所以实际上由于q˙(t)是一个pulse,先增后减的,我们在远场看到的求了导的结果就是先正后负的。
压强也一样。实际上三维的爆炸确实是这样的,外围先往外炸,冲击波过去之后会有负压。这部分低压空气立刻降温,里面的水分就凝结,形成condensation cloud。一些爆炸的照片上可以看到这种现象。
现在换一个情景。
假设有一个二维的圆形爆炸物,炸出来的东西只能在垂直于该炸弹的一维传播。那么,波动方程的解就是
ϕ(x,t)=2π∫0∞x2+s21f(t−cx2+s2)sds
令τ=t−x2+s2/c,把上式化为
2π∫0∞f(t−cx2+s2)dx2+s2
=2πc∫−∞t−x/cf(τ)dτ
=−2c∫−∞t−x/cq˙(τ)dτ
这个q是单位面积炸出的空间。可以看出来,这回,速度是v1=ϕ′(x,t)=21q˙(t−x/c),远场近场的行为是一样的。
现在,一个最有趣的东西是一维的爆炸物,它的冲击波可以在二维传播。这回,我们有
ϕ(x,t)=∫x∞r1f(t−cr)r2−x22rdr=2∫x∞f(t−cr)r2−x2dr
令τ=t−r/c,考虑f(t)=−q˙(t)/4π,r2−x2≈2x(r−x),我们有
ϕ(x,t)=2x2c∫−∞(t−x/c)f(τ)(ct−x)−cτdr=−2π1x2c∫−∞(t−x/c)q˙(τ)(t−x/c)−τdτ
对x求导,得到波速是这样的
v1(r,t)=2πc1x2c∫−∞(t−x/c)q¨(τ)(t−x/c)−τdτ
这东西是啥呢?实际上它可以看成是半阶导数。粗略地说,三维爆炸的远场行为是近场行为的一阶导,一维则是零阶导。半阶导怎么定义呢?可以用积分变换的方式:魔改一下拉普拉斯变换
(dtd)21F(t)= def π1∫−∞tt−τF˙(τ)dτ
所以v1实际上对应的是q˙的半阶导数。
远场行为如图。这玩意一直带一个不会消失的小尾巴,冲击波过了之后,相当于还有源源不断的向里的速度。这是非物理的,并且曾经在某些科学家的研究中造成过困难:一位科学家想要研究一个细长的裂缝导致的地震。他用了一维震源的近似,结果得到了这种非物理结果,害得他浪费大把时间在debug上。此事可引以为戒。
热传导方程
一、
∂t∂ϕ=κ∂x2∂2ϕ
可以用傅里叶变换来解:
ϕ(x,t)=∫−∞∞2πdkϕ~(k,t)eikx
∂t∂ϕ~=−κk2ϕ~
因此
ϕ(x,t)=∫−∞∞2πdkϕ~(k,t)eikx=∫−∞∞2πdkϕ~(k,0)eikx−κk2t
就得到一个原理上相当于格林函数的结果
ϕ(x,t)=∫−∞∞2πdk(∫−∞∞ϕ(ξ,0)eikξdξ)eikx−κk2t=∫−∞∞(∫−∞∞2πdkeik(x−ξ)−κk2t)ϕ(ξ,0)dξ=∫−∞∞G(x,ξ,t)ϕ(ξ,0)dξ
其中格林函数G(x,ξ,t)=∫−∞∞2πdkeik(x−ξ)−κk2t=4πκt1exp{−4κt1(x−ξ)2}就是所谓heat kernel。
二、
知乎上有一个民科,个人介绍是杜哈梅原理专家。我才知道杜哈梅就是duhamel,迪亚梅尔。我们先用一个例子看看这个原理的思路。
假设有一个一头受热(热源变温)的棍子,想解这个棍子的传热。这个方程非齐次,不好解。我们先考虑一个满足边界条件w(0,t)=1和w(x,0)=0的w,可以解出来是
w=θ(t){1−erf(2tx)}
其中erf是error function erf(x)=π2∫0xe−z2dz。
现在把每一个时刻的变温写成分立的
h(t)=∑nhnθ(t−tn)
或者
h(t)=h(0)+∫0th˙(τ)dτ=h(0)+∫0∞θ(t−τ)h˙(τ)dτ
由于传热方程是线性的,我们可以把每个时刻的解叠加起来,得到原非齐次问题的解
u(x,t)=∫0tw(x,t−τ)h˙(τ)dτ+h(0)w(x,t)=−∫0t(∂τ∂w(x,t−τ))h(τ)dτ=∫0t(∂t∂w(x,t−τ))h(τ)dτ
所以迪亚梅尔原理就是把非齐次方程变成一堆齐次的初值问题的叠加。
泊松和拉普拉斯方程
Schwinger’s trick
泊松方程的格林函数可以由g(r,r′)=∑nλn1φn(r)φn∗(r′)而写成
g(r,r′)=∫(2π)ndnkk2eik⋅(r−r′)
施温格的小技巧可以把这个积分算出来:
g(r,r′)=∫0∞ds∫(2π)ndnkeik⋅(r−r′)e−sk2=∫0∞ds(sπ)n(2π)n1e−4s1∣r−r′∣2=2nπn/21∫0∞dtt2n−2e−t∣r−r′∣2/4=2nπn/21Γ(2n−1)(4∣r−r′∣2)1−n/2=(n−2)Sn−11(∣r−r′∣1)n−2
注意到,对n=2的情况,这玩意发散。处理办法叫维度正规化,这是QFT里面的一种套路,其要义就是把火苗埋在地毯下面。在n=2附近展开
g(r,r′)=4π1(n/2−1)Γ(n/2)(1−(n/2−1)ln(41∣r−r′∣2)+O[(n−2)2])=4π1(n/2−11−2ln∣r−r′∣+⋯)
现在,我们只要把不取决于位置的1/(n/2−1)扫到毯子下面就行了。
惠更斯原理和基尔希霍夫近似
我们在中学都学过惠更斯的波的传播原理。这就是说,下一个时刻的波阵面,取决于上一个时刻波阵面上虚拟波源发出的次波的包络面。这是一个很朴素的想法,甚至21世纪的小学生也能想到。但是,波的传播还有另外一个特性:它能够相干叠加。惠更斯原理没有阐述这一点,而菲涅尔补全了这一部分,他的思路可以作为物理学家的思维方式的范本。
菲涅尔假设波前的任意一点Q在空间点P上贡献的复振幅为dψ(r,r′)。他根据以下直觉作出假设:
- 它正比于波阵面的面源;
- 它正比于r’处的次波源Q的复振幅;
- 它是球面波,dψ(r,r′)∝ReikR;
- 它是各向异性的。假设面元矢量和|r-r’|有夹角χ,那么dψ(r,r′)∝K(χ),即
ψ(r)=c∮Sψ(r′)K(χ)ReikRdS′
考虑波动方程的格林函数可以用波数作为参数写出来:
[−∇2−(c2ω2)]G(r,r′)=δ3(r−r′)
G±(r,r′)=4π∣r−r′∣1e±ik∣r−r′∣
假设在波源处只有向外的波。对于全空间,考虑到波动方程的微分算子是自伴的,使用Lagrange’s Identity(在这个特定的场景下,叫做格林公式):
∫Ω{G(r,r′)(∇r2+k2)ψ(r)−ψ(r)(∇r2+k2)G(r,r′)}dnx
=∫∂Ω{G(r,r′)∇rψ(r)−ψ(r)∇rG(r,r′)}⋅dSr
我们可以写
ψ(r′)=∫∂Ω{G(r,r′)(n⋅∇x)ψ(r)−ψ(r)∇rG(r,r′)}dSr,r′∈Ω
从这里其实可以看出来,要解波动方程,诺依曼和狄里赫雷边界条件只能给一个。我们假设从一个一维小孔(实际上就是小缝)里面射出平面波eikx(这意味着,真正产生光的点波源离小孔是很远的),
根据上图,其在小孔处的梯度是−ikeikx。在距离小孔的距离远大于波长的地方,可以近似
∇rG(r,r′)≈4πik∣r−r′∣2(r−r′)eik∣r−r′∣
代入,就知道
ψ(r′)≈4πik∫aperture ∣r−r′∣eik∣r−r′∣(1+cosθ)dSr
于是乎,因子K∝(1+cosθ)。这就是基尔希霍夫得出的修正因子。现在,新的“菲涅尔-基尔希霍夫”原理告诉我们两个惠更斯没有指出的事情:
- 波阵面往前产的波大于往回产的波。
- 分母的i给出一个90度的相位差,这是实验证实了的。对于这一点,如果不做计算,即使是惠更斯,恐怕也想不到吧。