下载文档原格式
第十章 偏微分方程数值解法偏微分方程问题，其求解十分困难。
除少数特殊情况外，绝大多数情况均难以求出精确解。
因此，近似解法就显得更为重要。
本章仅介绍求解各类典型偏微分方程定解问题的差分方法。
§1 差分方法的基本概念1.1 几类偏微分方程的定解问题椭圆型方程：其最典型、最简单的形式是泊松（Poisson ）方程),(2222y x f yu x u u =∂∂+∂∂=∆ 特别地，当0),(≡y x f 时，即为拉普拉斯（Laplace ）方程，又称为调和方程2222=∂∂+∂∂=∆yux u u Poisson 方程的第一边值问题为⎪⎩⎪⎨⎧Ω∂=Γ=Ω∈=∂∂+∂∂Γ∈),(),(),(),(),(2222y x y x u y x y x f y ux u y x ϕ 其中Ω为以Γ为边界的有界区域，Γ为分段光滑曲线，ΓΩ称为定解区域，),(y x f ，),(y x ϕ分别为Ω，Γ上的已知连续函数。
第二类和第三类边界条件可统一表示为),(),(y x u u y x ϕα=⎪⎪⎭⎫ ⎝⎛+∂∂Γ∈n 其中n 为边界Γ的外法线方向。
当0=α时为第二类边界条件， 0≠α时为第三类边界条件。
抛物型方程：其最简单的形式为一维热传导方程220(0)u ua a t x∂∂-=>∂∂ 方程可以有两种不同类型的定解问题：初值问题⎪⎩⎪⎨⎧+∞<<∞-=+∞<<-∞>=∂∂-∂∂x x x u x t x u a tu )()0,(,0022ϕ初边值问题221200,0(,0)()0(0,)(),(,)()0u ua t T x l t x u x x x lu t g t u l t g t t Tϕ⎧∂∂-=<<<<⎪∂∂⎪⎪=≤≤⎨⎪==≤≤⎪⎪⎩其中)(x ϕ，)(1t g ，)(2t g 为已知函数，且满足连接条件)0()(),0()0(21g l g ==ϕϕ边界条件)(),(),(),0(21t g t l u t g t u ==称为第一类边界条件。
第二十章 偏微分方程的数值解自然科学与工程技术中种种运动发展过程与平衡现象各自遵守一定的规律。
这些规律的定量表述一般地呈现为关于含有未知函数及其导数的方程。
我们将只含有未知多元函数及其偏导数的方程，称之为偏微分方程。
方程中出现的未知函数偏导数的最高阶数称为偏微分方程的阶。
如果方程中对于未知函数和它的所有偏导数都是线性的，这样的方程称为线性偏微分方程，否则称它为非线性偏微分方程。
初始条件和边界条件称为定解条件，未附加定解条件的偏微分方程称为泛定方程。
对于一个具体的问题，定解条件与泛定方程总是同时提出。
定解条件与泛定方程作为一个整体，称为定解问题。
§1 偏微分方程的定解问题各种物理性质的定常（即不随时间变化）过程，都可用椭圆型方程来描述。
其最典型、最简单的形式是泊松(Poisson)方程),(2222y x f y ux u u =∂∂+∂∂=∆ （1）特别地，当0),(≡y x f 时，即为拉普拉斯(Laplace)方程，又称为调和方程02222=∂∂+∂∂=∆yux u u （2）带有稳定热源或内部无热源的稳定温度场的温度分布，不可压缩流体的稳定无旋流动及静电场的电势等均满足这类方程。
Poisson 方程的第一边值问题为⎪⎩⎪⎨⎧Ω∂=Γ=Ω∈=∂∂+∂∂Γ∈),(|),(),(),(),(2222y x y x u y x y x f y uxu y x ϕ （3）其中Ω为以Γ为边界的有界区域，Γ为分段光滑曲线，ΓΩ 称为定解区域，),(),,(y x y x f ϕ分别为ΓΩ,上的已知连续函数。
第二类和第三类边界条件可统一表示成),(),(y x u n u y x ϕα=⎪⎭⎫⎝⎛+∂∂Γ∈ （4）其中n 为边界Γ的外法线方向。
当0=α时为第二类边界条件，0≠α时为第三类边界条件。
在研究热传导过程，气体扩散现象及电磁场的传播等随时间变化的非定常物理问题时，常常会遇到抛物型方程。
其最简单的形式为一维热传导方程)0(022>=∂∂-∂∂a xu a t u （5）方程（5）可以有两种不同类型的定解问题：初值问题（也称为Cauchy 问题）⎪⎩⎪⎨⎧+∞<<∞-=+∞<<∞->=∂∂-∂∂x x x u x t x u a t u )()0,(,0022ϕ （6）初边值问题⎪⎪⎪⎩⎪⎪⎪⎨⎧≤≤===<<<<=∂∂-∂∂Tt t g t l u t g t u x x u l x T t x ua t u 0),(),(),(),0()()0,(0,002122ϕ （7） 其中)(),(),(21t g t g x ϕ为已知函数，且满足连接条件)0()(),0()0(21g l g ==ϕϕ问题（7）中的边界条件)(),(),(),0(21t g t l u t g t u ==称为第一类边界条件。
第十六章 偏微分方程的数值解法科学研究和工程技术中的许多问题可建立偏微分方程的数学模型。
包含多个自变量的微分方程称为偏微分方程(partial differential equation)，简称PDE 。
偏微分方程问题，其求解是十分困难的。
除少数特殊情况外，绝大多数情况均难以求出精确解。
因此，近似解法就显得更为重要。
本章仅介绍求解各类典型偏微分方程定解问题的差分方法。
16.1 几类偏微分方程的定解问题一个偏微分方程的表示通常如下：(,,,,)x x x y y y x y A B C f x y Φ+Φ+Φ=ΦΦΦ (16.1.1) 式中，,,A B C 是常数，称为拟线性(quasilinear)数。
通常，存在3种拟线性方程： 双曲型(hyperbolic)方程：240B AC ->； 抛物线型(parabolic)方程：240B AC -=； 椭圆型(ellliptic)方程：240B AC -<。
16.1.2 双曲型方程最简单形式为一阶双曲型方程：0u ua t x∂∂+=∂∂ (16.1.2) 物理中常见的一维振动与波动问题可用二阶波动方程：22222u u a t x∂∂=∂∂ (16.1.3) 描述，它是双曲型方程的典型形式。
方程的初值问题为：2222200,(,0)()()t u uat x tx u x x u x x t ϕψ=⎧∂∂=>-∞<<+∞⎪∂∂⎪⎪=⎨⎪∂⎪=-∞<<+∞⎪∂⎩ (16.1.4)边界条件一般有三类，最简单的初边值问题为：2222212000,0(,0)(0,)(),(,)()0()t u ua t T x l t x u x lu t g t u l t g t t T ux x t ϕψ=⎧∂∂==<<<<⎪∂∂⎪⎪=≤⎪⎨==≤≤⎪⎪∂=-∞<<+∞⎪∂⎪⎩ (16.1.5)16.1.3 抛物型方程其最简单的形式为一维热传导方程：220(0)u ua a t x∂∂-=>∂∂ (16.1.8) 方程可以有两种不同类型的定解问题：(1) 初值问题：2200,(,0)()u ua t x t xu x x x ϕ⎧∂∂-=>-∞<<+∞⎪∂∂⎨⎪=-∞<<+∞⎩(16.1.6)(2) 初边值问题：221200,0(,0)()0(0,)(),(,)()0u ua t T x l t x u x x x l u t g t u l t g t t Tϕ⎧∂∂-=<<<<⎪∂∂⎪⎪=≤≤⎨⎪==≤≤⎪⎪⎩(16.1.7) 其中()x ϕ，1()g t ，2()g t 为已知函数，且满足连接条件：12(0)(0),()(0)g l g ϕϕ== (16.1.8)边界条件12(0,)(),(,)()u t g t u l t g t ==为第一类边界条件。
双重差分法（DID）安慰剂检验的做法：随机抽取500次？“安慰剂”（placebo）⼀词来⾃医学上的随机实验，⽐如要检验某种新药的疗效。
此时，可将参加实验的⼈群随机分为两组，其中⼀组为实验组，服⽤真药；⽽另⼀组为控制组，服⽤安慰剂（⽐如，⽆⽤的糖丸），并且不让参与者知道⾃⼰服⽤的究竟是真药还是安慰剂，以避免由于主观⼼理作⽤⽽影响实验效果。
双重差分法（DID）安慰剂检验的核⼼思想就是虚构处理组或者虚构政策时间进⾏估计，如果虚构情况下“伪政策虚拟变量”的系数依然显著，那么就说明原来的估计结果很有可能出现了偏误，我们的被解释变量y的变动很有可能是受到了其他政策或者随机性因素的影响。
说到虚构，那么⾃然是可以随机虚构，也可以不随机虚构（作者⾃⼰设定）。
当然，我更推荐的还是随机虚构处理组或者是政策时间的⽅法。
由于我们使⽤的数据基本都是“⼤N⼩T”型的短⾯板数据，所以随机虚构政策时间没什么意义，⽂献⼀般做法都是将政策年份统⼀提前2年或3年重新进⾏回归，看看政策虚拟变量系数是否依然显著。
我们更多地还是随机虚构处理组，具体做法就是随机选取个体作为处理组，重复500次或者1000次，看看“伪政策虚拟变量”的系数是否显著。
数据来源⽯⼤千等（2018）发表在《中国⼯业经济》的论⽂《智慧城市建设能否降低环境污染》使⽤DID⽅法评估了智慧城市建设对城市环境污染的影响，《中国⼯业经济》期刊官⽹公布了这篇论⽂使⽤的数据和代码。
接下来，我就使⽤这篇论⽂的数据，给⼤家分享⼀下双重差分法（DID）安慰剂检验中随机虚构处理组这种⽅法的Stata操作。
原⽂信息⽯⼤千,丁海,卫平,刘建江.智慧城市建设能否降低环境污染[J].中国⼯业经济,2018(06):117-135.随机虚构处理组的Stata操作双重差分法（DID）安慰剂检验的⼀般做法就是随机选取个体作为处理组，重复500次或者1000次，看看“伪政策虚拟变量”的系数是否显著。
在⽯⼤千等（2018）这篇论⽂中，处理组有32个城市，控制组有165个城市，所以我们需要从197个城市中随机选取32个城市作为“伪处理组”，假设这32个城市是智慧城市试点，其他城市为控制组，然后⽣成“伪政策虚拟变量”（交互项）进⾏回归。
精品文档引言主要讨论双曲性方程及双曲性方程组的差分方法。
从简单的一届线性双曲型方程开始，构造差分格式，分析其稳定性及其他性质，然后推广到一届线性双曲性方程组。
双曲方程与 椭圆方程，抛物方程的重要区别，是双曲方程具有特征和特征关系，其解对初值有局部依赖性质。
初值的函数性质（如间断，弱间断等）也沿特征传播，因而解一般无光滑性，迄今已发展许多逼近双曲方程的差分格式，这里只介绍常见的九种方法，讨论了各种求解方法，分析了其性质，最后对初边值问题及二维问题进行了讨论。
1 一阶线性常系数双曲型方程先考虑线性常系数方程[1]0=∂∂+∂∂xu a t u ,R x ∈,t>0 (1.1) 其中a 为给定常数，这是最简单的双曲型方程，一般称其为对流方程。
虽然（1.1）式非常简单，但是其差分格式的构造以及差分格式性质的讨论是讨论复杂的双曲型方程和方程组的基础。
它的差分格式可以推广到变系数方程，方程组以及拟线性方程和方程组。
对于方程（1.1）附以初始条件[1]u(x,0)=u 0(x), R x ∈ (1.2)在第一章中讨论了初值问题（1.1），（1.2）式的解，其解沿方程（1.1）的特征线[1]ε=-at x (1.3)是常数，并可表示为)()(),(00at x u u t x u -==ε下面讨论双曲性方程的应风格式，Lax-Friedrichs 格式，Lax-wendroff 格式，Courant-Friedrichs-Lewy 条件利用偏微分方程的特征线来构造有限差分格式，蛙跳格式，数值例子。
精品文档1.1 迎风格式迎风格式在实际计算中引起了普遍的重视，从而产生了很多好的方法和技巧。
迎风各式的 基本思想是简单的，就是在双曲型方程中关于空间偏导数用在特征线方向一侧的单边差商来代替，（1.1）式的迎风各式[1]是011=-+--+hu u au u nj n j n jn j τ, a>0 (1.4)的截断误差和稳定性:011=-+--+hu u au u nj n j n jn j τ, a>0+∂∂+∂∂+∂∂=-+3332221!31!21τττtu t u t u u un jn j/τ÷ +∂∂+∂∂+∂∂=-⇒+233221!31!21τττtu t u t u u u njn j ① n j n ju u 1-- +∂∂+∂∂-∂∂=333222!31!21h tu h x u h x u (两边乘于ha）,得 ⇒hu u anj n j 1+-=axu∂∂-!2a 22xu∂∂h +!3a 33xu ∂∂()420h h + ②①+②τn jn j u u -+1+hu u anj n j 1+-=tu∂∂+!2122t u ∂∂τ!3133tu ∂∂2τ+()30τ+a xu ∂∂-!2a 22xu∂∂h +!3a 33xu ∂∂()420h h +精品文档⎪⎭⎫ ⎝⎛∂∂+∂∂x u a t u+!2122t u∂∂τ-!2a 22xu∂∂h +!3a 33xu ∂∂2h +所以(,)j n T x t =!2122t u ∂∂τ-!2a 22xu∂∂h +截断误差为()h +τ0迎风格式对τ一阶精度,对h 一阶精度.当0,0h τ→→时(,)0j n T x t →,故迎风格式相容. 下面讨论迎风格式（1.4)的稳定性： 先把差分格式变化为便于计算的形式 n j n j u u -+1+()nj n j u u a 1--λ=0其中hτλ=网格式1+n j u =n j u -()n j n j u u a 1--λ令nj u =n u ikjh e则1+n v ikjh e =nv ikjh e -λa ()()h j ik n ikjh n e v e v 1--1+n v ikjh e =n v ikjh e -λa n v ikjh e +λa n v ikjh e ∙ikh e -1+n v=()ikh e a a -+-λλ1n v()k G ,τ=ikhae a -++λ1=()11-+-ikh e a λ=λa +1()1sin cos -+kh i kh =λa +1khcos -λa kh sin=λa +1()kh cos 1--λa i kh sin精品文档()2,k G τ=()[]2cos 11kh a --λ+kh a 222sin λ=()kh a cos 121--λ+22λa ()2cos 1kh -+kh a 222sin λ= ()kh a cos 121--λ+22λa +kh a 222sin λ=λa 221⨯-2cos 1kh-+22λa -22λa kh cos 2+22λa kh 2cos +22λa kh sin=λa 41-2sin 2kh+22λa -22λa kh cos 2+22λa =λa 41-2sin 2kh+222λa -222λa kh cos =λa 41-2sin 2kh+222λa ()kh cos 1- =λa 41-2sin 2kh +2⋅222λa 2cos 1kh-=λa 41-2sin 2kh+422λa 2sin2kh =λa 41-()λa +12sin 2kh当 1<λa 时原差分格式是稳定的。
