20190805个人笔记(川大)

官网课程视频地址

笔记主页

文献概要

主要讲解梯度流或守恒式方程
例如: Allen-Cahn: 或者 Cahn-Hilliard

Navier-stokes方程

相场模型(多相流)

模型问题

ϕt+Aϕ+N(ϕ)=0 \phi_t + A\phi + N(\phi)=0

方法一:需要解一串方程组
(ϕk)t+Aϕk+N(ϕ1,ϕ2)=0 (\phi_k)_t + A\phi_k + N(\phi_1, \phi_2,\cdots)=0
这种方法无法面对多维问题。

方法二:对时间离散

ϕn+1ϕnt+Aϕn+1+N(ϕn+1)=0 \frac{\phi^{n+1} - \phi^n}{\partial t}+ A\phi^{n+1} + N(\phi^{n+1})=0

模型来源

来源

ϕt=GδE[ϕ]δϕ \phi_t = -\mathcal G \frac{\delta E[\phi]}{\delta\phi}

根据链式法则
Et=(δE[ϕ]δϕ,ϕt) E_t = (\frac{\delta E[\phi]}{\delta\phi}, \phi_t )

G\mathcal G是一个正定算子,则是梯度流,此时
Et=(δE[ϕ]δϕ,GδE[ϕ]δϕ)0 E_t = - (\frac{\delta E[\phi]}{\delta\phi}, \mathcal G\frac{\delta E[\phi]}{\delta\phi} ) \leq 0

G\mathcal G反对称,则是Hamiltomain系统,此时
Et=0 E_t = 0

这里变分导数的定义

ddϵE[ϕ+ϵψ]ϵ=0=(δEδϕ,ψ) \frac{d}{d \epsilon} E[\phi + \epsilon \psi] |_{\epsilon=0}= (\frac{\delta E}{\delta \phi}, \psi)

E(ϕ)=E~(ϕ,ϕ) E(\phi) = \tilde E(\phi , \nabla \phi)
则有
Eϕ=E~ϕ+E~ϕ E_\phi = - \nabla \cdot\tilde E_\phi + \tilde E_{\nabla \phi}

例子1

根据G\mathcal G的不同选取,会得到不同的方程

G\mathcal G取正定算子
heat equation: E(ϕ)=Ω12ϕ2E(\phi)=\int_{\Omega} \frac{1}{2}|\nabla \phi|^{2} and G=I\mathcal{G}=I

Allen-Cahn: E(ϕ)=Ω(12ϕ2+14ϵ2(ϕ21)2)E(\phi)=\int_{\Omega}\left(\frac{1}{2}|\nabla \phi|^{2}+\frac{1}{4 \epsilon^{2}}\left(\phi^{2}-1\right)^{2}\right) and G=I\mathcal{G}=I

Cahn-Hilliard: E(ϕ)=Ω(ϵ2ϕ2+14ϵ(ϕ21)2)E(\phi)=\int_{\Omega}\left(\frac{\epsilon}{2}|\nabla \phi|^{2}+\frac{1}{4 \epsilon}\left(\phi^{2}-1\right)^{2}\right) and
G=Δ\mathcal{G}=-\Delta

Phase-field crystal: E(ϕ)=Ω(14ϕ4+α2ϕ2ϕ2+12Δϕ2)E(\phi)=\int_{\Omega}\left(\frac{1}{4} \phi^{4}+\frac{\alpha}{2} \phi^{2}-|\nabla \phi|^{2}+\frac{1}{2}|\Delta \phi|^{2}\right)
and G=Δ\mathcal{G}=-\Delta

L1L^{1} minimization: E(ϕ)=ΩϕE(\phi)=\int_{\Omega}|\nabla \phi| and G=I\mathcal{G}=I

G\mathcal G取反对称算子
Nonlinear Schrödinger equation:
E(ϕ)=Ω(12ϕ2+12F(ϕ2))E(\phi)=\int_{\Omega}\left(\frac{1}{2}|\nabla \phi|^{2}+\frac{1}{2} F\left(|\phi|^{2}\right)\right) and G=i\mathcal{G}=\mathrm{i}

一维情形
KDV equation: E(ϕ)=Ω(12xϕ2+ϕ3),G=xE(\phi)=\int_{\Omega}\left(\frac{1}{2}\left|\partial_{x} \phi\right|^{2}+\phi^{3}\right), \mathcal{G}=\partial_{x}

δEδϕ=ϕxxx+3ϕ2 \frac{\delta E}{\delta \phi} = - \phi_{xxx} + 3\phi^2

则有

ϕt=x(ϕxxx+3ϕ2)=xxxϕ+6ϕϕx \phi_t = - \partial_x(- \phi_{xxx} + 3\phi^2 ) = \partial_{xxx} \phi + 6\phi\phi_x

Allen-Cahn与Cahn-Hilliard能量稳定格式

引入一个新的变量μ\mu
{ϕt=Gμμ=Eϕ=Δϕ+F(ϕ) \begin{cases} {\phi_t} = -\mathcal G\mu\\ \mu =E_{\phi}= - \Delta \phi + F'(\phi) \end{cases}
边界条件ϕnΩ=μnΩ=0\frac{\partial \phi}{\partial \bm n}|_{\partial \Omega}=\frac{\partial \mu}{\partial \bm n}|_{\partial \Omega}= 0或者周期边界

Et=(Gμ,μ)0 E_t = -(\mathcal G \mu, \mu) \leq 0

证明数值格式的能量稳定性,原理是从原来格式中做相同的推算

全隐格式(Crank-Nicotson):

{ϕn+1ϕnΔt=Gμn+1+μn2μn+1+μn2=Δϕn+1+ϕn2+F(ϕn+1)F(ϕn)ϕn+1ϕn\begin{cases} \dfrac{\phi^{n+1} - \phi^n}{\Delta t} = -\mathcal G \dfrac{\mu^{n+1}+\mu^n}{2}\\[1.5em] \dfrac{\mu^{n+1}+\mu^n}{2} = -\Delta \dfrac{\phi^{n+1}+\phi^n}{2} +\dfrac{F(\phi^{n+1}) - F(\phi^n)}{\phi^{n+1} - \phi^n} \end{cases}

不难证出上式具有2阶无条件稳定,但是每次需要二阶非线性方程组,并且只能证明出当Δt\Delta t充分小的时候才能存在唯一解,该方法最早是杜强老师大约1991-1992年提出

凸分裂方法

对非线性项做分解
F(ϕ)=Fc(ϕ)Fe(ϕ) F(\phi) = F_c(\phi) - F_e(\phi)
并且FcF_cFeF_e都是凸函数,则可以构造下面的格式
{ϕn+1ϕnΔt=Gμn+1μn+1=Δϕn+1+Fc(ϕn+1)Fe(ϕn)\begin{cases} \dfrac{\phi^{n+1} - \phi^n}{\Delta t} = -\mathcal G \mu^{n+1}\\[1.5em] \mu^{n+1}= -\Delta \phi^{n+1} + F'_c(\phi^{n+1}) - F'_e(\phi^n) \end{cases}

注:一般而言FcF_c是能量耗散部分,FeF_e是能量递增格式

上式子分别与μn+1\mu^{n+1}Δt(ϕn+1ϕn)\Delta t(\phi^{n+1} - \phi^n)做内积

利用数学恒等式

2(ab,a)=(a,a)(b,b)+(ab,ab) 2(a-b, a) = (a, a) - (b,b) + (a-b, a-b)

FcF_cϕn+1\phi^{n+1}做泰勒展开

Fc(ϕn)=Fc(ϕn+1)+Fc(ϕn+1)(ϕnϕn+1)+12Fc(ξ)(ϕnϕn+1)2 F_c(\phi^n) = F_c(\phi^{n+1}) + F'_c(\phi^{n+1})(\phi^n - \phi^{n+1}) + \frac{1}{2}F''_c(\xi)(\phi^n - \phi^{n+1}) ^2
同理对FeF_eϕn\phi^{n}做泰勒展开

Fe(ϕn+1)=Fe(ϕn)Fe(ϕn+1)(ϕnϕn+1)+12Fe(η)(ϕnϕn+1)2 F_e(\phi^{n+1}) = F_e(\phi^{n}) - F'_e(\phi^{n+1})(\phi^n - \phi^{n+1}) + \frac{1}{2}F''_e(\eta)(\phi^n - \phi^{n+1}) ^2

两式子合并起来最后可以证明出该格式也是无条件稳定
缺点是: 只有一阶精度,也要解非线性方程组
优点是: 虽然也有非线性,但是非线性项是凸的,因此是存在唯一的,并且唯一解为凸泛函的能量最小解

证明: 只针对简单的情形G=I\mathcal G = I

ϕn+1ϕnΔt=Δϕn+1Fc(ϕn+1)+Fe(ϕn) \dfrac{\phi^{n+1} - \phi^n}{\Delta t}= \Delta \phi^{n+1} - F'_c(\phi^{n+1}) + F'_e(\phi^n)

定义Q(ϕ)=12Δtϕ2+12ϕ2+Fc(ϕ)gn(ϕ)Q(\phi) = \int\frac{1}{2\Delta t} |\phi|^2 + \frac{1}{2}|\nabla \phi|^2 + F_c(\phi) - g^n(\phi)

其中
gn(ϕ)=1Δtϕn+Fe(ϕn) g^n(\phi) = \frac{1}{\Delta t} \phi^n + F'_e(\phi^n)

凸泛函QQ的最小值就是变分导数为零的解

δQδϕϕ=ϕn+1=1Δtϕn+1Δϕn+1+Fc(ϕn+1)gn \frac{\delta Q}{\delta \phi}|_{\phi = \phi^{n+1}} = \frac{1}{\Delta t} \phi^{n+1} - \Delta \phi^{n+1} + F'_c(\phi^{n+1}) - g^n

该方法的特点

该方法最早提出来是1993年

最大缺点:

给一个二阶例子的构造F(ϕ)=(ϕ21)2F(\phi) = (\phi^2 -1)^2

Fc(ϕ)=14(ϕ4+1),Fe(ϕ)=12ϕ2F_c(\phi) = \frac{1}{4}(\phi^4+1), \quad F_e(\phi)=\frac{1}{2}\phi^2

则有Fc=ϕ3,Fe=ϕF'_c = \phi^3,\quad F'_e = \phi

{ϕn+1ϕnΔt=Gμn+1+μn2μn+1+μn2=Δϕn+1+ϕn2+A12(3ϕnϕn1)\begin{cases} \dfrac{\phi^{n+1} - \phi^n}{\Delta t} = -\mathcal G \dfrac{\mu^{n+1}+\mu^n}{2}\\[1.5em] \dfrac{\mu^{n+1}+\mu^n}{2} = -\Delta \dfrac{\phi^{n+1}+\phi^n}{2}+ A - \dfrac{1}{2}(3\phi^n - \phi^{n-1}) \end{cases}

其中
A=(ϕn+1)2(ϕn)22ϕn+1+ϕn2 A = \frac{(\phi^{n+1})^2 - (\phi^n)^2}{2} \frac{\phi^{n+1}+\phi^n}{2}

证明的过程中需要注意下面的拼凑

(3ϕnϕn1,ϕn+1ϕn)=(ϕn+1+ϕn(ϕn+12ϕn+ϕn1),ϕn+1ϕn)=(ϕn+1+ϕn(ϕn+1ϕn)+(ϕnϕn1)),ϕn+1ϕn)\begin{aligned} &(3\phi^n - \phi^{n-1}, \phi^{n+1} - \phi^n)\\ =&(\phi^{n+1}+\phi^n - (\phi^{n+1}-2\phi^n + \phi^{n-1}),\phi^{n+1}-\phi^n)\\ =&(\phi^{n+1}+\phi^n - (\phi^{n+1}-\phi^n) + (\phi^n-\phi^{n-1})),\phi^{n+1}-\phi^n) \end{aligned}

于是得到修正的能量是递减的,其中修正的能量定义为
E~(ϕn)=12ϕ2+14(ϕn+142ϕn+12)+14ϕn+1ϕn \tilde E(\phi^n) = \frac{1}{2}||\nabla \phi||^2 + \frac{1}{4}(||\phi^{n+1}||^4- 2||\phi^{n+1}||^2) + \frac{1}{4}|| \phi^{n+1}-\phi^n||

稳定化方法

从半隐格式出发
{ϕn+1ϕnΔt=Gμn+1μn+1=Δϕn+1+F(ϕn)\begin{cases} \dfrac{\phi^{n+1} - \phi^n}{\Delta t} = -\mathcal G \mu^{n+1}\\[1.5em] \mu^{n+1} = -\Delta \phi^{n+1} + F'(\phi^n) \end{cases}

优点: 计算快
缺点: 不稳定

改进办法: 加上稳定项

{ϕn+1ϕnΔt=Gμn+1μn+1=Δϕn+1+F(ϕn)\begin{cases} \dfrac{\phi^{n+1} - \phi^n}{\Delta t} = -\mathcal G \mu^{n+1}\\[1.5em] \mu^{n+1} = -\Delta \phi^{n+1} + F'(\phi^n) \end{cases}

能量修正
E(ϕ)=12ϕ2+Sϕ2+F(ϕ)Sϕ2dx E(\phi) = \int \frac{1}{2}|\nabla \phi|^2 + S|\phi|^2 + F(\phi) - S|\phi|^2 dx
或者
E(ϕ)=12(ϕ2+Sϕ2)(Sϕ2F(ϕ))dx E(\phi) = \int \frac{1}{2}|(\nabla \phi|^2 + S|\phi|^2 )-( S|\phi|^2-F(\phi))dx
Fe(ϕ)=Sϕ2F(ϕ)F_e(\phi) =S|\phi|^2 - F(\phi)
于是
Fe(ϕ)=2SF(ϕ) F''_e(\phi) = 2S - F''(\phi)
如果F(ϕ)L|F''(\phi)|\leq L,取SL2S\geq \frac{L}{2},则是特殊凸分裂格式,因此可以无条件稳定

早期能量是下面的格式
F(ϕ)=(1+ϕ)ln(1+ϕ)(1ϕ)ln(1ϕ)+θϕ2 F(\phi) = (1+\phi)\ln(1+\phi) - (1-\phi)\ln(1-\phi) + \theta \phi^2
物理上只关注ϕ(1,1)\phi\in(-1,1),

因此可以对能量修正,取
F~(ϕ)={14(ϕ21)2,ϕ1F(ϕ),ϕ>1. \tilde F(\phi)= \begin{cases} \frac{1}{4}(\phi^2 - 1)^2 , \quad & |\phi|\leq 1\\ F(\phi),\quad &| \phi| >1. \end{cases}

但是这种方法无法达到二阶

其实是现有了稳定化方法,再有凸分离方法,本质上这两种方法是有很大联系的

拉格朗日乘子法

F(ϕ)=14(ϕ21)2F(\phi) = \frac{1}{4}(\phi^2 - 1)^2

引入q=ϕ2q = \phi^2,则F(ϕ)=14q2F(\phi) = \frac{1}{4}q^2

{ϕt=Gμμ=Eϕ=Δϕ+qϕqt=2ϕϕt\begin{cases} \phi_t = -\mathcal G\mu \\ \mu = E_\phi = -\Delta \phi + q\phi\\ q_t = 2\phi\phi_t \end{cases}

数值格式为
{ϕn+1ϕn=ΔtGμn+1μn+1=Δϕn+1+qn+1ϕnqn+1qn=2ϕn(ϕn+1ϕn)\begin{cases} \phi^{n+1} - \phi^n =-\Delta t\mathcal G \mu^{n+1} \\ \mu^{n+1} = -\Delta \phi^{n+1} + q^{n+1}\phi^n\\ q^{n+1} - q^{n} = 2\phi^n(\phi^{n+1} - \phi^n) \end{cases}
上式分别与μn+1,ϕn+1ϕn,12qn+1\mu^{n+1}, \phi^{n+1}-\phi^n, \frac{1}{2}q^{n+1}做内积再联立可得修正的能量稳定,其中修正的能量定义为

E~(ϕn+1)=12ϕn+12+12qn+12 \tilde E({\phi^{n+1}}) = \frac{1}{2}|\nabla \phi^{n+1}|^2 + \frac{1}{2}||q^{n+1}||^2

优点:

结论:
做守恒型方程,用CN格式比较好,因为它没有多余的耗散项目
做耗散型方程,用BDF格式比较好

二阶格式简记:
对于CN格式:

ϕt\phi_t项,改写成
ϕn+1ϕnΔt \frac{\phi^{n+1}- \phi^n}{\Delta t}
对于ϕ\phi项,改写成
3ϕnϕn12 \frac{3\phi^{n} - \phi^{n-1}}{2}

对于BDF格式:
ϕt\phi_t项,改写成
3ϕn+14ϕn+ϕn12Δt \frac{3\phi^{n+1}- 4\phi^n+\phi^{n-1}}{2\Delta t}
对于ϕ\phi项,改写成
2ϕnϕn1 2\phi^{n} - \phi^{n-1}

不同的地方引入的辅助变量qq,CN格式12(qn+1+qn)\frac{1}{2}(q^{n+1}+q^n),BDF格式为qn+1q^{n+1}

注:拉格朗日乘子法是从方程出发的,直接对方程进行修正

傅里叶谱方法

{ϕΔϕ=fϕΩ=0,or ϕn=0,或周期边界 \begin{cases} \partial \phi - \Delta \phi = f\\ \phi|_{\partial \Omega}=0, \text{or } \frac{\partial \phi}{\partial n}| = 0, \text{或周期边界} \end{cases}
如果用傅里叶普方法
u=ukjeikxeijy u = \sum u_{kj}e^{ikx}e^{ijy}
于是
Δu=(k2+j2)ukjeikxeijy -\Delta u = \sum(k^2 + j^2)u_{kj}e^{ikx}e^{ijy}

作业

作业1

证明拉格朗日乘子法BDF-2格式的稳定性

作业2

写出Phase-field crystal的凸分裂格式格式

笔记主页