---
title: 6. Cahn-Hilliard equation
tags: Finite difference method
---

---

### Problem 5 (昆毅)

@Woodystick 1. 試著用 odesolver 解最簡單的 ode: $u_t = 3u, u(0)=1, t\in[0, 1]$ 試試. 我猜不管你 $\Delta t$ 取多少他幾乎都給你一樣的解. 2. 所以原題在 $t=0.1$ 時的解長怎樣? 能給我一個確定的嗎?
Ans:
1.嘗試過用固定步長的odesolver(MEBDF2)解最簡單的ODE有符合預期階數的的時間收斂性O(k^2), 但在這題上的時間收斂性還是不如預期, 所以我決定暫時放棄單獨討論時間收斂性
2.原題在 $t=0.1$ 時的解如下(原做法 k=h):
目測是可以收斂到真解的, 所以我下面就以原本的方法直接做了

@Woodystick (1) $e^{-0.1 \pi^4}\approx 5*10^{-5}$, 所以你 $h=0.04$ 算出來誤差 $10^{-4}$ 等於完全錯的 (2) 事實上最後一項 u_xx 也是線性的, 建議你第一步的 check 連那項一起做. (3) 時間上 O(k^4) 似乎完全無法 check? (4) 一般我們不會用 step, 直接用 t 就好. (5) true solution 是對的還錯的需要 check.

Cahn-Hilliard equation
$$
u_t = -0.1u_{xxxx}+(u^3-u)_{xx}, \quad x\in[-1, 1], \quad t\in[0, 1],
$$
with initial and boundary conditions
$$
u(x,0) = \cos(\pi x) - e^{-(2\pi x)^2}, \quad u(-1, t)=u(1, t) =-1, \quad u'(-1,t)=u'(1,t)=0.
$$

**Solution** (Code with Julia)
We use method of line to solve this problem.
But to check the accuracy and consistency of code,
we use the same method to solve a simple case first.

### Simple case:
$$
u_t = -0.1u_{xxxx}-u_{xx}, \quad x\in[-1, 1], \quad t\in[0, 0.1], \\
u(x,0) = \cos(\pi x) , \quad u(-1, t)=u(1, t) =- e^{-0.1\pi^4 t+\pi^2t}, \quad u'(-1,t)=u'(1,t)=0.
$$

with exact solution:
$$u(x,t)=e^{-0.1\pi^4 t+\pi^2t} \cos(\pi x) $$

We use central difference to approximate $u_{xxxx}$ and $u_{xx}$
$$
u_{xxxx} \approx
\frac{1}{h^4}(U_{i-2}-4U_{i-1}+6U_i-4U_{i+1}+U_{i+2}) \\
u_{xx} \approx
\frac{1}{h^2}(U_{i-1}-2U_i+1U_{i+1})
$$
With partition(ghost point):
$h=0.04, 0.02, 0.01$ (space step)
$N=\frac{2}{h}$ and $k=h$ (time step)
$x_i=ih-1$ with $i=-1 \sim N+1, (U_i=u(x_i,t))$
$t_i=ik$ with $i=0 \sim \frac{N}{2}$

$U_{0}=u(-1,t)=- e^{-0.1\pi^4 t+\pi^2t}= U_N=u(1,t),$
$u'(-1,t)=0 \Rightarrow \frac{U_{1}-U_{-1}}{2h}=0 \Rightarrow U_{-1}\approx U_{1},$
$u'(1,t)=0 \Rightarrow \frac{U_{N+1}-U_{N-1}}{2h}=0 \Rightarrow U_{N+1}\approx U_{N-1}$
(With truncation error $O(h^2)$)

then we get the system : $\frac{dU}{dt}= A\,U \text{, where } U=[U_1 \cdots U_{N-1}]^T$
$$\frac{dU}{dt}=\frac{1}{h^4} \begin{pmatrix}
-0.7+2h^2 & 0.4-h^2 & -0.1 & 0 & 0 & \cdots & 0 \\
0.4-h^2 & -0.6+2h^2 & 0.4-h^2 & -0.1 & 0 & \cdots & 0 \\
-0.1 & 0.4-h^2 & -0.6+2h^2 & 0.4-h^2 & -0.1 & \cdots & 0 \\
\vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \\
0 & \cdots &-0.1 & 0.4-h^2 & -0.6+2h^2 & 0.4-h^2&-0.1 \\
0 & \cdots &0&-0.1 & 0.4-h^2 & -0.6+2h^2 & 0.4-h^2\\
0 & \cdots & &0& -0.1 & 0.4-h^2 & -0.7+2h^2
\end{pmatrix}U
\\
\\
\text{with } U(x_i,0)=U^0_i=
\begin{pmatrix}
\cos (\pi x_1)\\
\vdots\\
\cos (\pi x_i)\\
\vdots\\
\cos (\pi x_N)
\end{pmatrix}
$$

Now we get an ode problem,
so we can use differential equation solvers of julia to solve it.
(In order to check the error of time step, we pick a method of fixed timestep. )
\
Method: MEBDF2 - The second order Modified Extended BDF method, which has improved stability properties over the standard BDF. Fixed timestep only.
(With truncation error $O(k^2)$)

h1=0.04
k1=h1
N1=Int64(round(2/h1, digits=0))-1
x1=(1:N1).*h1.-1
T=0.1
A_c = sparse([1;N1;2:N1-1;1:N1-1;2:N1;1:N1-2;3:N1], [1;N1;2:N1-1;2:N1;1:N1-1;3:N1;1:N1-2], 
    [-0.7+2*h1^2;-0.7+2*h1^2;(-0.6+2*h1^2)*ones(N1-2);(0.4-h1^2)*ones(N1-1);
        (0.4-h1^2)*ones(N1-1);(-0.1)*ones(N1-2);(-0.1)*ones(N1-2)])

u0 = cos.(pi.*x1)
M=t->[-0.4*exp(-0.1*pi^4*t+pi^2*t)+h1^2*exp(-0.1*pi^4*t+pi^2*t);0.1*exp(-0.1*pi^4*t+pi^2*t);zeros(N1-4);
    0.1*exp(-0.1*pi^4*t+pi^2*t);-0.4*exp(-0.1*pi^4*t+pi^2*t)+h1^2*exp(-0.1*pi^4*t+pi^2*t)]
f(u,p,t) = (A_c*u+p(t))/h1^4
tspan = (0.0,T)
prob = ODEProblem(f,u0,tspan,M)
sol = solve(prob,MEBDF2(),dt=k1)

we can get the solution:

### error analysis:
because we know the exact solution of this problem is
$$u(x,t)=e^{-0.1\pi^4 t+\pi^2t} \cos(\pi x) $$
we can get the error $||U(x_i,0.1)-u(x_i,0.1)||$
which should be $O(h^2)+O(k^2)$ at $t=1$

change k=h :
| | k=h=0.04 | k=h=0.02 | k=h=0.01 |
| -------- | -------- | -------- |-------- |
|1 norm|0.000788282|0.000197293| 4.93196e-5|
|2 norm|0.000743397 | 0.00018607 | 4.65372e-5|
|$\infty$ norm|0.000944258|0.000236418|5.91609e-5|

$$ R_\infty=\frac{||U2-U3||_\infty}{||U1-U2||_\infty}=0.25042023212843845 \approx \frac{1}{4} $$
(R值結果符合預期)

不同步長的解於T=0.1時的圖(幾乎重疊):
上圖在中間部分放大的結果:
總上, 隨著空間步長的遞減, 解有收斂到真實解的趨勢, 此方法應該可以接受, 故以相同方法做本題實際要求的情形

### True case:
Now consider
$$
u_t = -0.1u_{xxxx}+(u^3-u)_{xx}, \quad x\in[-1, 1], \quad t\in[0, 1],
$$
with initial and boundary conditions
$$
u(x,0) = \cos(\pi x) - e^{-(2\pi x)^2}, \quad u(-1, t)=u(1, t) =-1, \quad u'(-1,t)=u'(1,t)=0.
$$

We use central difference to approximate $u_{xxxx}$ and $u_{xx}$
\
$$
u_{xxxx} \approx
\frac{1}{h^4}(U_{i-2}-4U_{i-1}+6U_i-4U_{i+1}+U_{i+2})
$$
$$
u_{xx} \approx
\frac{1}{h^2}(U_{i-1}-2U_i+1U_{i+1})
$$
With partition(ghost point):
$h=0.04, 0.02, 0.01$ (space step)
$N=\frac{2}{h}$ and $k=h$ (time step)
$x_i=ih-1$ with $i=-1 \sim N+1 , (U_i=u(x_i,t))$
$t_i=ik$ with $i=0 \sim \frac{N}{2}$

$U_{0}=u(-1,t)=- 1= U_N=u(1,t),$
$u'(-1,t)=0 \Rightarrow \frac{U_{1}-U_{-1}}{2h}=0 \Rightarrow U_{-1}\approx U_{1},$
$u'(1,t)=0 \Rightarrow \frac{U_{N+1}-U_{N-1}}{2h}=0 \Rightarrow U_{N+1}\approx U_{N-1}$
(With truncation error $O(h^2)$)

then we get the system : $\frac{dU}{dt}= A\,U+B\,U^3+p(t) \text{, where } U=[U_1 \cdots U_{N-1}]^T , \, U^3=[U_1^3 \cdots U_{N-1}^3]^T$

$$\frac{dU}{dt}=\frac{1}{h^4} \begin{pmatrix}
-0.7+2h^2 & 0.4-h^2 & -0.1 & 0 & 0 & \cdots & 0 \\
0.4-h^2 & -0.6+2h^2 & 0.4-h^2 & -0.1 & 0 & \cdots & 0 \\
-0.1 & 0.4-h^2 & -0.6+2h^2 & 0.4-h^2 & -0.1 & \cdots & 0 \\
\vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \\
0 & \cdots &-0.1 & 0.4-h^2 & -0.6+2h^2 & 0.4-h^2&-0.1 \\
0 & \cdots &0&-0.1 & 0.4-h^2 & -0.6+2h^2 & 0.4-h^2\\
0 & \cdots & &0& -0.1 & 0.4-h^2 & -0.7+2h^2
\end{pmatrix}U
$$
$$
+\frac{1}{h^2} \begin{pmatrix}
-2 & 1 & 0 & \cdots & 0 \\
1 & -2 & 1 & \cdots & 0 \\
\vdots & \ddots & \ddots & \ddots & \vdots \\
0 & \cdots & 1 & -2 & 1 \\
0 & \cdots & 0 & 1 & -2
\end{pmatrix}U^3
+\frac{1}{h^4}
\begin{pmatrix}
-0.4\\
0.1\\
0\\
\vdots\\
0\\
0.1\\
-0.4
\end{pmatrix}
\\
\\
\text{with } U(x_i,0)=U^0_i=
\begin{pmatrix}
\cos (\pi x_1)-e^{-4 \pi^2x_1^2}\\
\vdots\\
\cos (\pi x_i)-e^{-4 \pi^2x_i^2}\\
\vdots\\
\cos (\pi x_{N-1})-e^{-4 \pi^2x_{N-1}^2}
\end{pmatrix}
$$

and like we do above, we use MEBDF2 of julia to solve it.
(With truncation error $O(k^2)$)

T=1 ##t \in [0,1]
h1=0.04
k1=h1
N1=Int64(round(2/h1, digits=0))-1
x1=(1:N1).*h1.-1
A_c = sparse([1;N1;2:N1-1;1:N1-1;2:N1;1:N1-2;3:N1], [1;N1;2:N1-1;2:N1;1:N1-1;3:N1;1:N1-2], 
    [-0.7+2*h1^2;-0.7+2*h1^2;(-0.6+2*h1^2)*ones(N1-2);(0.4-h1^2)*ones(N1-1);
        (0.4-h1^2)*ones(N1-1);(-0.1)*ones(N1-2);(-0.1)*ones(N1-2)])
B_c = sparse([1;N1;2:N1-1;1:N1-1;2:N1], [1;N1;2:N1-1;2:N1;1:N1-1], 
    [2;2;(2)*ones(N1-2);(-1)*ones(N1-1);
        (-1)*ones(N1-1)])
f(u,p,t)=(A_c*u-h1^2*B_c*u.^3+p(t))/h1^4
u0 = cos.(pi.*x1).-exp.(-4*pi^2*x1.^2)
M=t->[-0.4;0.1;zeros(N1-4);0.1;-0.4]
tspan = (0.0,T)
prob = ODEProblem(f,u0,tspan,M)
sol = solve(prob,MEBDF2(),dt=k1)

解隨時間變化圖:

不同步長的解於T=0.1時的圖(幾乎重疊):
上圖在中間部分放大的結果(有收斂趨勢):

不同步長的解於T=1時的圖(也有收斂趨勢, 肉眼可看出不同是因為y軸其實被放大了):
所以本題實際解(T=1時)應很接近上圖情形

---

### Problem 6 (浩恩)

Cahn-Hilliard equation
$$
u_t = -0.02u_{xxxx}+(u^3-u)_{xx}, \quad x\in[-4, 4], \quad t\in[0, 100],
$$
with initial and boundary conditions
$$
u(x,0) = -1 + 10\left(e^{-35(x+2)^2} + e^{-11 x^2} + e^{-7 (x-2)^2}\right), \quad u(-4, t)=u(4, t) =-1, \quad u'(-4,t)=u'(4,t)=0.
$$

---
**Solution**

#### Simple case
$$
u_t = -u_{xxxx}+u_{xx}, \quad x\in[-4, 4], \quad t\in[0, 1],\\
u(x,0) = \cos(\frac{\pi}{2} x), \quad u(-4, t)=u(4, t) =e^{-(\frac{\pi^4}{16}+\frac{\pi^2}{4})t}, \quad u'(-4,t)=u'(4,t)=0.
$$
with exact solution $u(x,t)=cos(\frac{\pi}{2} x)e^{-(\frac{\pi^4}{16}+\frac{\pi^2}{4})t}$

We divide the domain $x=[-4-h=x_{-1},-4=x_0,-4+h=x_1,\cdots,4-h=x_{N-1},4=x_N,4+h=x_{N+1}]$ where $N=\frac{8}{h}.$

From the B.C., we have $U^n_0=U^n_N=e^{-(\frac{\pi^4}{16}+\frac{\pi^2}{4})nk},U^n_{-1}=U^n_1,U^n_{N+1}=U^n_{N-1}$ for $n\leq\frac{1}{k}.$ ($k$ is the time step size)

Using Crank-Nicolson and central difference method, we have the discrete form :

$\frac{U^{n+1}_j-U^n_j}{k} = \frac{-1}{2h^4}[(U^{n+1}_{j+2}-4U^{n+1}_{j+1}+6U^{n+1}_j-4U^{n+1}_{j-1}+U^{n+1}_{j-2})+(U^n_{j+2}-4U^n_{j+1}+6U^n_j-4U^n_{j-1}+U^n_{j-2})]\\
\quad\quad\quad\quad+\frac{1}{2h^2}[(U^{n+1}_{j+1}-2U^{n+1}_j+U^{n+1}_{j-1})+(U^n_{j+1}-2U^n_j+U^n_{j-1})]$

We can get the following system :
$(h^4I+A)U^{n+1}=(h^4I-A)U^n+C,$ where

$A=\begin{pmatrix}
7a-2b & -4a+b & a & 0 & 0 & \cdots & 0 \\
-4a+b & 6a-2b & -4a+b & a & 0 & \cdots & 0 \\
a & -4a+b & 6a-2b & -4a+b & a & \cdots & 0 \\
\vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \\
0 & \cdots &a & -4a+b & 6a-2b & -4a+b&a \\
0 & \cdots &0&a & -4a+b & 6a-2b & -4a+b\\
0 & \cdots & &0& a & -4a+b & 7a-2b
\end{pmatrix}$ and $C=(U_0^{n+1}+U_0^n)\begin{pmatrix}
4a-b\\
-a\\
0\\
\vdots\\
0\\
-a\\
4a-b
\end{pmatrix}$ with $a = \frac{k}{2}$ and $b=\frac{-kh^2}{2}.$

Here is the exact solution with $T=0,1$

Error :
Since we use Crank-Nicolson and central difference method, the order of $h$ and $k$ should both be 2.

|$(h,k)$|$\|u-$exact$\|_\infty$|
|---|---|
|$(\frac{1}{25},\frac{1}{50})$|6.165415326915636e-5|
|$(\frac{1}{50},\frac{1}{100})$|1.5420919750746008e-5|
|$(\frac{1}{100},\frac{1}{200})$|3.855638904410263e-6|
|$(\frac{1}{200},\frac{1}{400})$|9.611951430971361e-7|
From the table above, we can check that the truncation error is $O(h^2)+O(k^2).$

using LinearAlgebra
using Plots
using LaTeXStrings
using SparseArrays

# discrete space
h = 1/50;
N = Int64(8/h);
x = [-4+h:h:4-h;];

# time step
k = 1/100; # dt
T = 1;  # end time
n = Int64(ceil(T/k)); # total time step number

# initial u(x,0) and solution of t=1,2,...,T
U = zeros(N-1,T+1);
U[:,1] = [cos(pi/2*x[i]) for i in 1:length(x)]; # u(x,0)

# construct matrix
a = k/2;
b = -k*(h^2)/2;
A = sparse([1;2:N-2;N-1;1:N-2;2:N-1;1:N-3;3:N-1;],[1;2:N-2;N-1;2:N-1;1:N-2;3:N-1;1:N-3;],
            [7*a-2*b;(6*a-2*b)*ones(N-3);7*a-2*b;(-4*a+b)*ones(N-2);(-4*a+b)*ones(N-2);a*ones(N-3);a*ones(N-3);]);
C = sparsevec([1;2;N-2;N-1;],[4*a-b;-a;-a;4*a-b;]);

# iteration
u = U[:,1];
for i = 1:n
    bd1 = exp((-pi^4/16-pi^2/4)*i*k); # boundary value at t = ik
    bd0 = exp((-pi^4/16-pi^2/4)*(i-1)*k); # boundary value at t = (i-1)k
    u = ((h^4)*I+A)\Array(((h^4)*I-A)*u+(bd0+bd1)*C);

    if i % Int64(1/k) == 0
        U[:,Int64(i*k)+1] = Array(u);
    end
end

# error
norm(U[:,end]-exact,Inf)

#### True case
Now get back to our original problem :
$$
u_t = -0.02u_{xxxx}+(u^3-u)_{xx}, \quad x\in[-4, 4], \quad t\in[0, 100],
$$
with initial and boundary conditions
$$
u(x,0) = -1 + 10\left(e^{-35(x+2)^2} + e^{-11 x^2} + e^{-7 (x-2)^2}\right), \quad u(-4, t)=u(4, t) =-1, \quad u'(-4,t)=u'(4,t)=0.
$$

Using the same method above, we have the discrete form :

$\frac{U^{n+1}_j-U^n_j}{k} = \frac{-0.02}{2h^4}[(U^{n+1}_{j+2}-4U^{n+1}_{j+1}+6U^{n+1}_j-4U^{n+1}_{j-1}+U^{n+1}_{j-2})+(U^n_{j+2}-4U^n_{j+1}+6U^n_j-4U^n_{j-1}+U^n_{j-2})]\\
\quad\quad\quad\quad-\frac{1}{2h^2}[(U^{n+1}_{j+1}-2U^{n+1}_j+U^{n+1}_{j-1})+(U^n_{j+1}-2U^n_j+U^n_{j-1})]\\
\quad\quad\quad\quad+\frac{1}{2h^2}\{[(U^{n+1}_{j+1})^3-2(U^{n+1}_j)^3+(U^{n+1}_{j-1})^3]+[(U^n_{j+1})^3-2(U^n_j)^3+(U^n_{j-1})^3]\}$

We can get the following system :
$(h^4I+A)U^{n+1}-B(U^{n+1})^3=(h^4I-A)U^n+B(U^n)^3+C,$ where

$A=\begin{pmatrix}
7a-2b & -4a+b & a & 0 & 0 & \cdots & 0 \\
-4a+b & 6a-2b & -4a+b & a & 0 & \cdots & 0 \\
a & -4a+b & 6a-2b & -4a+b & a & \cdots & 0 \\
\vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \\
0 & \cdots &a & -4a+b & 6a-2b & -4a+b&a \\
0 & \cdots &0&a & -4a+b & 6a-2b & -4a+b\\
0 & \cdots & &0& a & -4a+b & 7a-2b
\end{pmatrix},\\
B=\begin{pmatrix}
-2b & b & 0 & \cdots & 0 \\
b & -2b & b & \cdots & 0 \\
\vdots & \ddots & \ddots & \ddots & \vdots \\
0 & \cdots & b & -2b & b \\
0 & \cdots & 0 & b & -2b
\end{pmatrix},\quad C=\begin{pmatrix}
-8a\\
2a\\
0\\
\vdots\\
0\\
2a\\
-8a
\end{pmatrix}$ and $a=\frac{0.02k}{2}, b=\frac{kh^2}{2}.$

To solve the nonlinear system, we use Newton method where

$J(U^k)=h^4I+A+\begin{pmatrix}
6b(U^k_1)^2 & -3b(U^k_2)^2 & 0 & \cdots & 0 \\
-3b(U^k_1)^2 & 6b(U^k_2)^2 & -3b(U^k_3)^2 & \cdots & 0 \\
\vdots & \ddots & \ddots & \ddots & \vdots \\
0 & \cdots & -3b(U^k_{N-3})^2 & 6b(U^k_{N-2})^2 & -3b(U^k_{N-1})^2 \\
0 & \cdots & 0 & -3b(U^k_{N-2})^2 & 6b(U^k_{N-1})^2
\end{pmatrix}$

## setting
# discrete space
h = 1/100;
N = Int64(8/h);
x = [-4+h:h:4-h;];

# time step
k = 1/1200; # dt
T = 100;  # end time
n = Int64(ceil(T/k)); # total time step number

#safetime = [Int64(1/k)*[0.1;0.5;1;10;25;50;75;100];];
safetime = [10;60;600;Int64(1/k)*[1;10;25;50;75;100];]

# initial u(x,0) and solution of t = safetime
U = zeros(N-1,length(safetime)+1);
U[:,1] = [-1+10*(exp(-35*((x[i]+2)^2))+exp(-11*(x[i]^2))+exp(-7*((x[i]-2)^2))) for i in 1:length(x)]; # u(x,0)

## construct matrix
a = 0.02*k/2;
b = k*(h^2)/2;

A = sparse([1;2:N-2;N-1;1:N-2;2:N-1;1:N-3;3:N-1;],[1;2:N-2;N-1;2:N-1;1:N-2;3:N-1;1:N-3;],
            [7*a-2*b;(6*a-2*b)*ones(N-3);7*a-2*b;(-4*a+b)*ones(N-2);(-4*a+b)*ones(N-2);a*ones(N-3);a*ones(N-3);]);
B = sparse([1:N-1;1:N-2;2:N-1;],[1:N-1;2:N-1;1:N-2;],[-2*b*ones(N-1);b*ones(N-2);b*ones(N-2);]);
C = sparsevec([1;2;N-2;N-1;],[-8*a;2*a;2*a;-8*a;]);

## solve the nonlinear system for each time step
u = U[:,1]; # U^0
tol = 1e-10;
for i =1:n
    err = 1;
    RHS = (h^4*I-A)*u + B*(u.^3) + C;
    u = ones(N-1);
    while err > tol
        F = (h^4*I+A)*u - B*(u.^3) - RHS;
        J = h^4*I + A + sparse([1:N-1;1:N-2;2:N-1;],[1:N-1;2:N-1;1:N-2;],[6*b*u.^2;-3*b*u[2:end].^2;-3*b*u[1:end-1].^2;]);
        temp = J\Array(F);
        u = u - temp;
        err = norm(temp,Inf);
    end
    if i in safetime
        U[:,findfirst(isequal(i), safetime)+1] = Array(u);
    end
end

We only choose some time point to save : $t=0,\frac{1}{120},\frac{1}{20},\frac{1}{2},1,10,25,50,75,100$

Set $h=\frac{1}{100}$ and $k=\frac{1}{1200}$, the numerical result :

The result with different $(h,k)$ at $T=100$ :

They are almost the same, it seems that our result is the answer we want.