---
title: 7. Kuramoto-Sivashinsky equation
tags: Finite difference method
---

---

### Problem 7 (銘瑜)

Kuramoto-Sivashinsky equation
$$
u_t = -0.08u_{xxxx}-u_{xx}+u u_{x}, \quad x\in[-2, 2], \quad t\in[0, 20],
$$
with initial and periodic boundary conditions
$$
u(x,0) = 1 + 0.5e^{-40 x^2}.
$$

---

$h=\frac{4}{N-1}$ is grid size in space.
$k$ is grid size in time.
Backward Euler(中間的點):
$$
\frac{U^{n+1}_i-U^n_i}{k}=-0.08\frac{U^{n+1}_{i+2}-4U^{n+1}_{i+1}+6U^{n+1}_{i}-4U^{n+1}_{i-1}-U^{n+1}_{i-2}}{h^4}-
\frac{U^{n+1}_{i+1}-2U^{n+1}_{i}+U^{n+1}_{i-1}}{h^2}
+U^{n+1}_{i}\frac{U^{n+1}_{i+1}+U^{n+1}_{i-1}}{2h}
$$

Equals to
$$
U^{n}_{i}=0.08k\frac{U^{n+1}_{i+2}-4U^{n+1}_{i+1}+6U^{n+1}_{i}-4U^{n+1}_{i-1}-U^{n+1}_{i-2}}{h^4}
+k\frac{U^{n+1}_{i+1}-2U^{n+1}_{i}+U^{n+1}_{i-1}}{h^2}
-U^{n+1}_{i}k\frac{U^{n+1}_{i+1}+U^{n+1}_{i-1}}{2h}+U^{n+1}_{i}
$$

Periodic boundary conditions(最前面2個點跟最後面兩個點):

$$
U^{n}_{1}=0.08k\frac{U^{n+1}_{3}-4U^{n+1}_{2}+6U^{n+1}_{1}-4U^{n+1}_{N}-U^{n+1}_{N-1}}{h^4}
+k\frac{U^{n+1}_{2}-2U^{n+1}_{1}+U^{n+1}_{N}}{h^2}
-U^{n+1}_{1}k\frac{U^{n+1}_{2}+U^{n+1}_{N}}{2h}+U^{n+1}_{1}
$$

$$
U^{n}_{2}=0.08k\frac{U^{n+1}_{4}-4U^{n+1}_{3}+6U^{n+1}_{2}-4U^{n+1}_{1}-U^{n+1}_{N}}{h^4}
+k\frac{U^{n+1}_{3}-2U^{n+1}_{2}+U^{n+1}_{1}}{h^2}
-U^{n+1}_{2}k\frac{U^{n+1}_{3}+U^{n+1}_{1}}{2h}+U^{n+1}_{2}
$$

$$
U^{n}_{N-1}=0.08k\frac{U^{n+1}_{1}-4U^{n+1}_{N}+6U^{n+1}_{N-1}-4U^{n+1}_{N-2}-U^{n+1}_{N-3}}{h^4}
+k\frac{U^{n+1}_{N}-2U^{n+1}_{N-1}+U^{n+1}_{N-2}}{h^2}
-U^{n+1}_{N}k\frac{U^{n+1}_{N}+U^{n+1}_{N-2}}{2h}+U^{n+1}_{N-1}
$$

$$
U^{n}_{N}=0.08k\frac{U^{n+1}_{2}-4U^{n+1}_{1}+6U^{n+1}_{N}-4U^{n+1}_{N-1}-U^{n+1}_{N-2}}{h^4}
+k\frac{U^{n+1}_{1}-2U^{n+1}_{N}+U^{n+1}_{N-1}}{h^2}
-U^{n+1}_{N}k\frac{U^{n+1}_{1}+U^{n+1}_{N-1}}{2h}+U^{n+1}_{N}
$$

To solve nonlinear system as above,use Newton method. Therefore,input initial condition $u_0$ we will get the solution of the next time step $u_1$,what we need to do is keeping this process $20/k$ times.

#### ERROR

##### Local truncation error

$$
\tau=u_t+0.08u_{xxxx}+u_{xx}+u u_{x}-\\(\dfrac{U_{i}^{n+1}-U_{i}^{n}}{k}+0.08(\dfrac{U_{i-2}^{n+1}-4U_{i-2}^{n+1}+6U_{i}^{n+1}-4U_{i+1}^{n+1}+U_{i+2}^{n+1}}{h^4})+\\\dfrac{U_{i+1}^{n+1}-2U_{i}^{n+1}+U_{i-1}^{n+1}}{h^2}-U_{i}^{n+1}\dfrac{(U_{i+1}^{n+1})-(U_{i-1}^{n+1})}{2h})=O(k+h^2)
$$

Here we compare the two-norm between several size of h and k at $t=0.016$.
let $h_1=4/100$ , $h_2 = 4/200$ , $h_3 = 4/400$ , $h_4 = 4/800$ $h_5 = 4/1600$ and $k_1 = h_1^2$, $k_2 = h_2^2$, $k_3 = h_3^2$ , $k_4 = h_4^2$ , $k_5 = h_5^2$

$$
R_{1} = \frac{\lVert U_{h_1}^{k_1} - U_{h_2}^{k_2} \rVert_2}{\lVert U_{h_2}^{k_2} - U_{h_3}^{k_3} \rVert_2}=\frac{0.0025856968353127688}{0.0006961738397084752}=3.714153976822135
$$

$$
R_{2} = \frac{\lVert U_{h_2}^{k_2} - U_{h_3}^{k_3} \rVert_2}{\lVert U_{h_3}^{k_3} - U_{h_4}^{k_4} \rVert_2}=\frac{0.0006961738397084752}{0.0001771030195276322}=3.9308976298953264 \approx 2^2
$$

$$
R_3 = \frac{\lVert U_{h_3}^{k_3} - U_{h_4}^{k_4} \rVert_2}{\lVert U_{h_4}^{k_4} - U_{h_5}^{k_5} \rVert_2}=\frac{0.0001771030195276322}{4.446117216283617e-05} =3.9833187231097695$ \approx 2^2
$$

It seems converge.

#### Code

import numpy as np
import numpy.linalg as npLA
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve
import matplotlib.pyplot as plt
import matplotlib.animation as animation

n=100
k = 0.001
iteraion=20000

x=np.arange(-2,2,4/n)
h=x[1]-x[0]
t=np.arange(iteraion+1)*k
u=np.zeros( ( n, iteraion+1 ) )
u_0=1+0.5*np.exp(-40*x**2)#initialcondition

alpha = 0.08
diagonals = (k/h**4)*np.array([h**2-4*alpha , alpha , alpha , h**2-(4*alpha) , 6*alpha-2*h**2+1/(k/h**4) , h**2-4*alpha , alpha , alpha ,  h**2-4*alpha ])
LA = diags( diagonals , [-n+1 , -n+2 , -2,  -1, 0 , 1, 2, n-2, n-1] , shape=(n,n) , format='csc')

u0=u_0

#Newton method objective function
def Newton(n,k,h,u0):
    u_n=u0.reshape(n, 1)
    for i in range(10):
        u0=u0.reshape(n, 1)
        diagonals = (k/2/h)*np.array([1,1,1,1])
        F=diags( diagonals,[-n+1,n-1,-1,1] , shape=(n,n) , format='csc')
        object_F=LA*u0-F.dot(u0.reshape(n, 1))*u0.reshape(n, 1)-u_n
        u0=u0.reshape(1,n)
        #Jocobian
        diagonals=(k/2/h)*np.array([-u0,np.roll(u0,1),-u0 , np.roll(u0,-1)-np.roll(u0,1) , u0 ])
        dF = diags(diagonals , [-n+1,n-1,-1,0, 1] , shape=(n,n) , format='csc')
        J = LA-dF
        u0=u0-spsolve(J,object_F)
        if(npLA.norm(object_F))<1e-8:
            break
    return u0
fig = plt.figure()
ims = []
im = plt.plot(x,u)
ims.append(im)
for i in range(iteraion):
    u0=Newton(n,k,h,u0)
    if 0==i%100: 
        im = plt.plot(x,u0,'b')
        ims.append(im)

#### Result:

This is an example where h=4/100,k=0.001,
There are 200 pictures in this gif,and the gap between then is 0.1sec,which means from start to end it will take exactly 20secs.

藍線為$h=4/100$, $k=0.001$,橘線為$h=4/200$,$k=0.00025$,,綠線為 $h=4/400$,$k=0.0000625$

T=1

T=3

T=5

T=7

T=9

T=20

可以看出前面幾秒是收斂的，但超過9秒，已經開始出現不小的差異，
順帶一題在$h=4/50$, $k=0.002$下，大約在7秒就已經和上面的橘藍綠線出現差異
，而要畫出T=20有收斂的結果必須把網格點切得很細。

#### 2020/5/8
藍線為$h=4/400$, $k=0.001$,橘線為$h=4/800$,$k=0.001$,綠線為 $h=4/1600$,$k=0.001$
,紅線為$h=4/400$,$k=0.0001$

藍橘綠三條線重疊，由此推斷之前T=20看起來不收斂並非空間上的點切的不夠細。

藍線為$h=4/100$, $k=0.00001$,橘線為$h=4/100$,$k=0.000005$,,綠線為 $h=4/100$,$k=0.0000025$
將時間網格點切得非常細，得到看起來收斂的結果。

### Problem 8 (尚文)

Kuramoto-Sivashinsky equation
$$
u_t = -0.07u_{xxxx}-u_{xx}+u u_{x}, \quad x\in[-2, 2], \quad t\in[0, 20],
$$
with initial and periodic boundary conditions
$$
u(x,0) = 1 + 0.5e^{-40 x^2}.
$$

#### Solution

(Note: $\alpha=0.07$)

I tried Backward Euler and Crank-Nicolson with Newton's method and failed. Therefore, I used linearization that does linear computation.

Grid design:
$$
x_i = -2 + ih, x_0 = -2, x_N= 2-h, h=4/N
$$

The periodic boundary condition is enforced by identifying $U_{-1}= U_N, U_{-2} = U_{N-1}, U_{N+1}= U_0, U_{N+2} = U_1$. This can be implemented in `python` using `np.roll`, which allows the grid variables being indexed by $\frac{\mathbb{Z}}{(N+1)\mathbb{Z}}$ instead of $\mathbb{Z}$.

#### 3-step "linearized" backward Euler

Using standard BE, the non-linear relation is
$$
\frac{U^{n+1}i-U^{n}_i}{k}=-\alpha\frac{U{i-2}^{n+1}-4U_{i-1}^{n+1}+6U_i^{n+1}-4U_{i+1}^{n+1}+U_{i+2}^{n+1}}{h^4}-\frac{U^{n+1}{i-1}-2U_i^{n+1}+U^{n+1}_{i+1}}{h^2}+U^{n+1}_i\frac{U^{n+1}{i+1}-U^{n+1}_{i-1}}{2h}
$$

I substituded $U^{n+1}_i$ in the non-linear part with an extrapolation in $t$ -
$$
U^{n+1}_i\approx 2 U_i^n- U^{n-1}_i
$$

Thus for each $n\ge 2$, we solve
$$
F_i(\vec{U^{n-1}},\vec{U^n},\vec{U^{n+1}})=U^{n+1}_i-U^{n}_i+k\alpha\frac{U_{i-2}^{n+1}-4U_{i-1}^{n+1}+6U_i^{n+1}-4U_{i+1}^{n+1}+U_{i+2}^{n+1}}{h^4}+k\frac{U^{n+1}_{i-1}-2U_i^{n+1}+U^{n+1}_{i+1}}{h^2}-k(\frac12 U^{n}_i-\frac32 U_i^{n-1})\frac{U^{n+1}_{i+1}-U^{n+1}_{i-1}}{2h}=0
$$

We can turn this system of equation into the matrix form
$$
A(\vec{U^{n-1}},\vec{U^n})\vec{U^{n+1}}=\vec{U^n}
$$

And for $n=1$, I used FE to explicitly calculate
$$
U^1_i=U_i^0-k\alpha\frac{U_{i-2}^{0}-4U_{i-1}^{0}+6U_i^{0}-4U_{i+1}^{0}+U_{i+2}^{0}}{h^4}
$$

The code is as follows

def Matrix_3steps(alpha, h,k,U_past, U_pastpast):
    leg = np.shape(U_past)[0]
    N = leg - 1
    h = 4/N
    
    #CIRCULANT TERM
    c_2 = k*alpha/h**4
    c_1 = -k*4*alpha/h**4 + k/h**2
    c_0 = 1 + k*6*alpha/h**4 - 2*k/h**2

    J = peanut.diags([c_1, c_2, c_2, c_1, c_0, c_1, c_2, c_2, c_1], [N, N-1, 2, 1, 0, -1, -2, -N+1, -N], format = 'csr')
    
    #NON-CIRCULANT TERM
    Non = peanut.csr_matrix( (leg, leg) )
    
    ro = np.arange(leg)
    ro_1 = np.roll(ro,-1)
    rop1 = np.roll(ro,1)
    
    for i in range(leg):
        fake_unp1 = 2*U_past[i,0] - 1*U_pastpast[i,0]
        #superdiagonal
        Non[i,rop1[i]] = -fake_unp1*(k/(2*h))
        #subdiagonal
        Non[i,ro_1[i]] = fake_unp1*(k/(2*h))
        
    return J + Non

def KS_3steps_20secs(alpha, grid, k): # to play 20 seconds, we must iterate for 20/k times, and that k must be s.t. 20/k is integer
    gridnum = np.shape(grid)[0]
    h = grid[1] - grid[0]
    
    #t=0k
    init_condit = []
    for i in range(gridnum):
        init_condit.append(init_condit_gen(grid[i]))
    init_condit = np.reshape(np.array(init_condit), (gridnum,1))
    lamb_u = [init_condit]
    
    #t=1k
    ro = np.arange(gridnum)
    ro_2 = np.roll(ro,-2)
    ro_1 = np.roll(ro,-1)
    rop1 = np.roll(ro,1)
    rop2 = np.roll(ro,2)
    
    
    init_next = np.zeros((gridnum,1))
    for i in range(gridnum):
        U_i_2 = init_condit[ro_2[i],0]
        U_i_1 = init_condit[ro_1[i],0]
        U_i = init_condit[i,0]
        U_ip1 = init_condit[rop1[i],0]
        U_ip2 = init_condit[rop2[i],0]
        
        init_next[i,0] = U_i - (k/h**4)*0.07*(U_i_2-4*U_i_1+6*U_i-4*U_ip1+U_ip2) - (k/h**2)*(U_i_1-2*U_i+U_ip1) + (k/(2*h))*(U_ip1-U_i_1)
    lamb_u.append(init_next)
    
    #t=nk
    n_t = int(20/k) 
    for i in tqdm(range(1, n_t - 2)): #扣掉 t=0k, t=1k
        M = Matrix_3steps(alpha, h, k, lamb_u[i], lamb_u[i-1])
        b = nutella.spsolve( M, lamb_u[i])
        lamb_u.append( columnize(b) )
    
    return lamb_u

##### Local truncation error

$u_t\approx\frac{U^{n+1}_i-U^{n}_i}{k}$
the error of this term is $O(k)$

$u_{xxxx}\approx\frac{U_{i-2}^{n+1}-4U_{i-1}^{n+1}+6U_i^{n+1}-4U_{i+1}^{n+1}+U_{i+2}^{n+1}}{h^4}$
the error of this term is $O(h^2)$

$u_{xx}\approx\frac{U^{n+1}_{i-1}-2U_i^{n+1}+U^{n+1}_{i+1}}{h^2}$
the error of this term is $O(h^2)$

$uu_x\approx U^{n+1}_i\frac{U^{n+1}_{i+1}-U^{n+1}_{i-1}}{2h}$
the error of this term is $O(h^2)$

<<May 18>>

So the expected error is $O(k+h^2)$. In particular, if we set $k=\gamma h^2$, then the error is expected to be $O(h^2)$

Computing $R=\frac{\|U_{h}-U_{h/2}\|}{\|U_{h/2}-U_{h/4}\|}$, we expect the number to be $2^2$
At $t=0.1$,
`R1 = 4.402289261921088`
`R2 = 4.990612405846505`
`Rinf = 5.933569038913571`

@AAAaEAfsRJSoei9DGIWSfg For the results shown above, what is $\Delta t$ and $\Delta h$ that you choose? Also, try to refine the mesh to see if the values are getting closer to $4$.

<< May 25 >>
I tried 5 kinds of grids and multiple t:
5 kinds of grids -
(1) `h=4/128, k=1/64`
(2) `h=4/256, k=1/256`
(3) `h=4/512, k=1/1024`
(4) `h=4/1024, k=1/4096`
(5) `h=4/2048, k=1/16384`

Thus there are three kinds of R's
$R^1=\frac{\|U_{h=128}-U_{h=256}\|}{\|U_{h=256}-U_{h=512}\|}$ (128-256-512)

$R^2=\frac{\|U_{h=256}-U_{h=512}\|}{\|U_{h=512}-U_{h=1024}\|}$ (256-512-1024)

$R^3=\frac{\|U_{h=512}-U_{h=1024}\|}{\|U_{h=1024}-U_{h=2048}\|}$ (512-1024-2048)

Where the indices indicate the grid size, not the norm.

For 2-grid norm, the result is as follows:

|t=|1/16|2/16|3/16|4/16|5/16|6/16|7/16|8/16|9/16|
|---|---|---|---|---|---|---|---|---|---|---|
|$R_2^1$|1125.7816| 57.7525| 19.6439| 19.3507| 20.0505| 20.9901| 22.4576| 24.5522| 27.3145|
|$R_2^2$|3.7949| 3.3521| 4.6897| 5.2423| 4.9906| 4.7162| 4.4964| 4.2511| 3.9116|
|$R_2^3$|3.3966| 3.7541| 3.8665| 3.406| 3.1142| 2.9859| 2.886| 2.774| 2.6459|

where the R associated with the coarsest grids ($R_2^1$) have extremely large numbers, but $R_2^2$ fluctuates around $2^2$, while $R_2^3$ is even smaller and is below $2^2$.

##### Result

I tried three grids
(1) `h=4/256, k=1/256`
(2) `h=4/512, k=1/1024`
(3) `h=4/1024, k=1/4096`

(https://imgur.com/a/w7Bx2ix)

For $L^2$ norm,
20 seconds -

5 seconds -

---

### Problem 9 (沛聖)

Kuramoto-Sivashinsky equation
$$
u_t = -0.06u_{xxxx}-u_{xx}+u u_{x}, \quad x\in[-2, 2], \quad t\in[0, 20],
$$
with initial and periodic boundary conditions
$$
u(x,0) = 1 + 0.5e^{-40 x^2}.
$$

#### Numerical Algorithm

Given $N\in \mathbf{N}$. Consider the grid points $(x_i,t_n) = ( -2+(i-1)h , nk),\,i=1,\cdots,N$, where $h=\frac{4}{N}$ (mesh width), $k=\Delta t$ (time step).

Let $U_i^n \approx u(x_i,t_n)$ denote the numerical approximation at grid point $(x_i,t_n)$.

Note that the periodic boundary condition implies that $U_{k+N}^n = U_{k}^n,\, \forall n\geq 0,\, \forall k\in \mathbf{N}$

Use Backward Euler method, we obtain
$$
D_t U_i^{n+1} = -\alpha D_{xxxx} U_i^{n+1} - D_{xx} U_i^{n+1} + U_i^{n+1} D_x U_i^{n+1}
$$
where
1. $\alpha=0.06$
1. $D_t$ denotes the $1^{st}$-order backward difference w.r.t. $t$ ,
1. $D_{xxxx}$ denotes the $4^{th}$-order central difference w.r.t. $x$.
1. $D_{xx}$ denotes the $2^{nd}$-order central difference w.r.t. $x$.
1. $D_{x}$ denotes the $1^{st}$-order central difference w.r.t. $x$.

\begin{align*}
\implies \frac{U_i^{n+1}-U_i^n}{k} &= -\alpha\left(\frac{U_{i-2}^{n+1} - 4U_{i-1}^{n+1} + 6U_i^{n+1} -4U_{i+1}^{n+1} + U_{i+2}^{n+1} }{h^4}\right)
- \left(\frac{U_{i-1}^{n+1} - 2U_i^{n+1} + U_{i+1}^{n+1} }{h^2}\right)
+ U_i^{n+1} \left(\frac{U_{i+1}^{n+1}-U_{i-1}^{n+1} }{2h}\right)\\
&= -\left( \frac{ \alpha U_{i-2}^{n+1}+(h^2-4\alpha)U_{i-1}^{n+1} + (6\alpha-2h^2)U_i^{n+1} +(h^2-4\alpha)U_{i+1}^{n+1} +\alpha U_{i+2}^{n+1} }{ h^4} \right)+ U_i^{n+1} \left(\frac{U_{i+1}^{n+1}-U_{i-1}^{n+1} }{2h}\right)
\end{align*}

Let $r = \frac{k}{h^4},\, s=\frac{k}{2h}$, then this can be written as
\begin{align*}
& r\left( \alpha U_{i-2}^{n+1}+(h^2-4\alpha)U_{i-1}^{n+1} + (6\alpha-2h^2+\frac{1}{r})U_i^{n+1} +(h^2-4\alpha)U_{i+1}^{n+1} +\alpha U_{i+2}^{n+1} \right) - s U_i^{n+1}\left(U_{i+1}^{n+1}-U_{i-1}^{n+1}\right) = U_i^n\\
\iff & AU^{n+1} - F(U^{n+1}) = U^n
\end{align*}

where
$$
A=r\begin{pmatrix}
6\alpha-2h^2+1/r & h^2-4\alpha & \alpha & & && \alpha & h^2-4\alpha\\
h^2-4\alpha & 6\alpha-2h^2+1/r & h^2-4\alpha & \alpha & && & \alpha\\
\alpha & h^2-4\alpha & 6\alpha-2h^2+1/r & h^2-4\alpha & \alpha && &\\
&\ddots&\ddots&\ddots&\ddots&\ddots&&\\
& &\ddots&\ddots&\ddots&\ddots&\ddots&\\
&&& \alpha & h^2-4\alpha & 6\alpha-2h^2+1/r & h^2-4\alpha &\alpha\\
\alpha &&&& \alpha & h^2-4\alpha & 6\alpha-2h^2+1/r & h^2-4\alpha\\
h^2-4\alpha & \alpha &&&& \alpha & h^2-4\alpha & 6\alpha-2h^2+1/r
\end{pmatrix}
$$
and
$$
F(U) =
\begin{pmatrix}
F_1(U) \\ F_2(U) \\ \vdots \\ F_{N-1}(U)\\F_N(U)
\end{pmatrix}
= s
\begin{pmatrix}
U_1(U_2-U_N) \\ U_2(U_3-U_1) \\ \cdots \\ U_{N-1} (U_N-U_{N-2}) \\ U_N( U_1 - U_{N-1})
\end{pmatrix}
$$

Let $G(U) = AU - F(U)$, then $U^{n+1}$ can be found by solving
$$
G^{(n)}(U) := G(U) - U^n = \vec{0}
$$
with Newton's method, given that $U^n$ is known.

In order to run Newton's method, we need to compute Jacobian of $G^{(n)}(U) = J_{G^{(n)}}(U)$.

In particular,
\begin{align*}
J_{G^{(n)}} (U) &= J_G(U) = \left[ \frac{\partial G_i}{\partial U_{(j)}} \right]_{i,j=1}^N,\,\quad U_{(j)}=\text{j-th coordinate of vector } U\\
&= A - \left[ \frac{\partial F_i}{\partial U_{(j)}} \right]_{i,j=1}^N\\
&= A - s
\begin{pmatrix}
U_2 - U_N & -U_2 & & & & U_N\\
U_1 & U_3 - U_1 & -U_3 & & &\\
& U_2 & U_4 - U_2 & -U_4 & &\\
& & \ddots & \ddots & \ddots & \\
& & & \ddots & \ddots &\ddots\\
-U_1 & & & & U_{N-1} & U_1 - U_{N-1}
\end{pmatrix}
\end{align*}

##### Algorithm

1. $X = \{-2+(i-1)h | i=1,\dots,N,\, h=4/N \},\quad k=\Delta t$
1. $U^0 = 1+0.5 \exp (-40X^2)$
1. $\text{for }i=1,\dots,N \\
U^k \leftarrow U^{k-1}\\
\text{while not convergent:}\\
U^k \leftarrow U^k - J_{G^{(k-1)}}(U^k)^{-1} G^{(k-1)}(U^k)\\
\text{endwhile}\\
\text{endfor}$

##### Remark
\begin{align*}
J_{G^{(n)}} (U) &= J_G(U) = \left[ \frac{\partial G_i}{\partial U_{(j)}} \right]_{i,j=1}^N,\,\quad U_{(j)}=\text{j-th coordinate of vector } U\\
&= A - \left[ \frac{\partial F_i}{\partial U_{(j)}} \right]_{i,j=1}^N\\
&= A - s
\begin{pmatrix}
U_2 - U_N & -U_2 & & & & U_N\\
U_1 & U_3 - U_1 & -U_3 & & &\\
& U_2 & U_4 - U_2 & -U_4 & &\\
& & \ddots & \ddots & \ddots & \\
& & & \ddots & \ddots &\ddots\\
-U_1 & & & & U_{N-1} & U_1 - U_{N-1}
\end{pmatrix}
\end{align*}

##### Code

import numpy as np
import numpy.linalg as npLA
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve
from scipy.optimize import newton
import matplotlib.pyplot as plt

#### Input ####
# L: absolute value of right/left endpoints. The spacial domain is [-L,L].
# N: uniformly divide spacial domain [-L,L] into N subintervals.
# dt: time step
# T : maximal time 

#### Output ####
# x: mesh grid
# u: solution. u[i,n] = u(x_i , t_n) = u(i*dx,n*dt)


def KSsolve(L,N,dt,T):
    h = 2*L/N
    k = dt
    
    # iter_max: maximal iteration of time t. iter_max*k < T
    iter_max = int(T/k)

    x = np.arange(-L,L,h)
    t = np.arange(iter_max+1)*k
    u = np.zeros( ( N, iter_max+1 ) )

    u[:,0] = 1+0.5*np.exp(-40*x**2)

    r = k/h**4
    s = k/2/h
    alpha = 0.06
    diagonals = r*np.array([h**2-4*alpha , alpha , alpha , h**2-(4*alpha) , 6*alpha-2*h**2+1/r , h**2-4*alpha , alpha , alpha ,  h**2-4*alpha ])
    A = diags( diagonals , [-N+1 , -N+2 , -2,  -1, 0 , 1, 2, N-2, N-1] , shape=(N,N) , format='csc')
    
    # Newton's method
    for n in range(iter_max):
        u0 = u[:,n]
        u_old = u0
        for i in range(100):
            update = spsolve(df(u_old,A,s),func(u_old,u0,A,s))
            u_new = u_old - update
            u_old  = u_new
            if npLA.norm(func(u_new,u0,A,s)) < 1e-8:
                break
        print("n=" , n , "/",iter_max ," i=",i)
        u[:,n+1] = u_new
        
    return x,t,u

def func(u,u0,A,s):
    f = s*np.multiply( u, np.roll(u,-1)-np.roll(u,1) )
    y = A*u -f - u0
    return y

def df(u,A,s):
    N = A.shape[0]
    diagonals = s*np.array([-u , u , np.roll(u,-1)-np.roll(u,1) , -np.roll(u,-1) , np.roll(u,1)])
    dF = diags(diagonals , [-N+1 , -1 , 0 , 1 , N-1] , shape=(N,N) , format='csc')
    Jac = A-dF
    return Jac

#### Grid convergence study
##### Local Truncation Error
Suppose that $u(x,t)$ is the true solution of our K-S equation, and let $u_i^n = u(x_i,t_n)$,

then the local truncation error is
\begin{align*}
&\tau(x_i,t_n) \\
=& D_t u_i^n + \alpha D_{xxxx} u_i^{n+1} +D_{xx} u_i^{n+1} - u_i^{n+1} D_x u_i^{n+1}\\
=& (\frac{1}{2} u_{tt} + u_{xxt}+ u_{xxxxt} - u_x u_t - u u_{xt})k + (\frac{1}{12} u_{xxxx} + \frac{1}{6} u_{xxxxxx} - \frac{1}{6}uu_{xxx})h^2 +o(k+h^2),\quad \text{at } (x,t)=(x_i,t_n)
\end{align*}

##### Estimating Error
Let $U_h^k$ denote the numerical solution with grid size $h$ and time step $k$.

Choose $h_1=4/100$ , $h_2 = 4/200$ , $h_3 = 4/400$ and $k_1 = h_1^2$, $k_2 = h_2^2$, $k_3 = h_3^2$ and $t=5$,

then
$$
R = \frac{\lVert U_{h_1}^{k_1} - U_{h_2}^{k_2} \rVert_2}{\lVert U_{h_2}^{k_2} - U_{h_3}^{k_3} \rVert_2}=3.9693473788076306 \approx 2^2
$$

##### Remark
If we choose $k=h^2$, then the local truncation error $$\tau(x_i,t_n) = Ck+Dh^2 = (C+D)h^2$$

##### Long time behavior of L2 energy
Assuming that the domain is $[-L,L]$, it is shown, numerically, that the energy of long-time solutions satisfies
$$
\lim\sup_{t\to\infty} \lVert u(\cdot,t) \rVert_{L_2} = \lim\sup_{t\to\infty} \left[ \int_{-L}^L u(x,t)^2\,dx \right]^{1/2} = O(L^{1/2})
$$

Set $h=2L/100,\,\Delta t=0.005$. Consider different cases of problem with parameter $L=1,2,3,4,5$, we compute the energy, $\lVert u(\cdot,t) \rVert_{L_2}$, of the numerical solution $u(x,t)$ and plot the energy functions with respect to different time span.

1. $0\leq t \leq 10$
1. $0 \leq t \leq 100$
Observe that after $t\approx 15$, solution $u$ for the case $L=5$ become periodic in space.


1. $0 \leq t \leq 1000$

###### Ref:
1. [Kuramoto-Sivashinsky equation. Encyclopedia of Mathematics.](https://www.encyclopediaofmath.org/index.php/Kuramoto-Sivashinsky_equation)
1. [An in-depth numerical study of the two-dimensional Kuramoto–Sivashinsky equation](https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.2014.0932)

#### Numerical Result
Here I choose $N=100$, time step $\Delta t = 0.01$ and maximal time $T=20$.

L = 2
N = 100
dt = 0.01
T = 20

(x,t,u) = KSsolve(L,N,dt,T)


---

### Problem 10 (暐軒)

Kuramoto-Sivashinsky equation
$$
u_t = -0.05u_{xxxx}-u_{xx}+u u_{x}, \quad x\in[-2, 2], \quad t\in[0, 20],
$$
with initial and periodic boundary conditions
$$
u(x,0) = 1 + 0.5e^{-40 x^2}.
$$

---

Answer:
## Algorithm
* Partition:
$x_1=-2, \ x_{N+1}=2,\ x_i=-2+ih$ with $h=\dfrac{4}{N},\ i=1,...,N$
k is the time step , $t_n=nk$
\
$U_i^n=u(x_i,t_n)$ for $i=1,...,N$

* Discretization:
$\dfrac{U_{i}^{n+1}-U_{i}^{n}}{k}=-0.05\dfrac{U_{i-2}^{n+1}-4U_{i-2}^{n+1}+6U_{i}^{n+1}-4U_{i+1}^{n+1}+U_{i+2}^{n+1}}{h^4}-\dfrac{U_{i+1}^{n+1}-2U_{i}^{n+1}+U_{i-1}^{n+1}}{h^2}+U_{i}^{n+1}\dfrac{(U_{i+1}^{n+1})-(U_{i-1}^{n+1})}{2h}$
Let F($\bf{U^{n+1}}$)=
$\dfrac{U_{i}^{n+1}-U_{i}^{n}}{k}+0.05(\dfrac{U_{i-2}^{n+1}-4U_{i-2}^{n+1}+6U_{i}^{n+1}-4U_{i+1}^{n+1}+U_{i+2}^{n+1}}{h^4})+\dfrac{U_{i+1}^{n+1}-2U_{i}^{n+1}+U_{i-1}^{n+1}}{h^2}-U_{i}^{n+1}\dfrac{(U_{i+1}^{n+1})-(U_{i-1}^{n+1})}{2h}$
\
Find the jacobian matrix of F
\
$J(\bf{U})=A-g(U)$
\
linear term A:
\begin{pmatrix}
1+6\dfrac{0.05k}{h^4}-\dfrac{2k}{h^2}&\dfrac{k}{h^2}-4\dfrac{0.05k}{h^4}&\dfrac{0.05k}{h^4}&...&0&\dfrac{0.05k}{h^4}&\dfrac{k}{h^2}-4\dfrac{0.05k}{h^4}\\
\dfrac{k}{h^2}-4\dfrac{0.05k}{h^4} & 1+6\dfrac{0.05k}{h^4}-\dfrac{2k}{h^2} & \dfrac{k}{h^2}-4\dfrac{0.05k}{h^4}&\dfrac{0.05k}{h^4}&...&0&\dfrac{0.05k}{h^4}\\
\dfrac{0.05k}{h^4}&\dfrac{k}{h^2}-4\dfrac{0.05k}{h^4}&1+6\dfrac{0.05k}{h^4}-\dfrac{2k}{h^2}&\dfrac{k}{h^2}-4\dfrac{0.05k}{h^4}&\dfrac{0.05k}{h^4}&\dots&0\\
&&...\\
&&...\\
&&...\\
&&...\\
\dfrac{0.05k}{h^4}&0&...&\dfrac{0.05k}{h^4}&\dfrac{k}{h^2}-4\dfrac{0.05k}{h^4}&1+6\dfrac{0.05k}{h^4}-\dfrac{2k}{h^2}&\dfrac{k}{h^2}-4\dfrac{0.05k}{h^4}\\
\dfrac{k}{h^2}-4\dfrac{0.05k}{h^4}&\dfrac{0.05k}{h^4}&0&...&\dfrac{0.05k}{h^4}&\dfrac{k}{h^2}-4\dfrac{0.05k}{h^4}&1+6\dfrac{0.05k}{h^4}-\dfrac{2k}{h^2}
&
\end{pmatrix}
\
nonlinear term $g(\bf{U})$
\
\begin{pmatrix}
k\dfrac{U_{2}^{n+1}-U_{N}^{n+1}}{2h}&\dfrac{kU_{1}^{n+1}}{2h}&...&...&...&...&\dfrac{-kU_{1}^{n+1}}{2h}\\
\dfrac{-kU_{2}^{n+1}}{2h}&k\dfrac{U_{3}^{n+1}-U_{1}^{n+1}}{2h}&\dfrac{kU_{2}^{n+1}}{2h}&0&...&...&0\\
0&\dfrac{-kU_{3}^{n+1}}{2h}&k\dfrac{U_{4}^{n+1}-U_{2}^{n+1}}{2h}&\dfrac{kU_{3}^{n+1}}{2h}&...&...&0\\
&&...\\
&&...\\
&&...\\
&&...\\
&&...\\
\dfrac{kU_{N}^{n+1}}{2h}&...&...&...&...&\dfrac{kU_{N}^{n+1}}{2h}&k\dfrac{U_{1}^{n+1}-U_{N-1}^{n+1}}{2h}\\
\end{pmatrix}
\
Apply the Netown's method
\
${\bf{U^{n+1}}}=\bf{U^{n}}-{J(U^n)}^{-1}F(\bf{U^{n}})$

## Convergence
* Local truncation error
$u_t+0.05u_{xxxx}+u_{xx}+u u_{x}-(\dfrac{U_{i}^{n+1}-U_{i}^{n}}{k}+0.05(\dfrac{U_{i-2}^{n+1}-4U_{i-2}^{n+1}+6U_{i}^{n+1}-4U_{i+1}^{n+1}+U_{i+2}^{n+1}}{h^4})+\dfrac{U_{i+1}^{n+1}-2U_{i}^{n+1}+U_{i-1}^{n+1}}{h^2}-U_{i}^{n+1}\dfrac{(U_{i+1}^{n+1})-(U_{i-1}^{n+1})}{2h})\\=O(k)+O(h^2)$

* The order of error in time is O($k$).
為確定時間網格上的收斂，下面會固定總時間T，取不同time step
觀察求出的解圖形上是否重合及 R值是否為2

$$
R=\frac{||U_k-U_\frac{k}{2}||}{||U_\frac{k}{2}-U_\frac{k}{4}||}=2^p,
$$
where p is the order of error

The total number of points :N=100
Total time T=10.
we have 3 diferent time step
k1=0.0001 ; k2=0.00005 ;k3=0.000025
這裡有三條對應不同time step的線，藍線、橘線和綠線
\
\
$R=\frac{||U_k-U_\frac{k}{2}||}{||U_\frac{k}{2}-U_\frac{k}{4}||}=2.004002932524179$
\
Total time T=15.
we have 3 diferent time step
k1=0.0001 ; k2=0.00005 ;k3=0.000025
\
\
$R=\frac{||U_k-U_\frac{k}{2}||}{||U_\frac{k}{2}-U_\frac{k}{4}||}=1.880826364596718$
\
若將時間的網格切得更細
k1=0.00005; k2=0.000025; k3=0.0000125
\
\
$R=\frac{||U_k-U_\frac{k}{2}||}{||U_\frac{k}{2}-U_\frac{k}{4}||}=1.9390715785188408$
\
Total time T=20.
we have 3 diferent time step
k1=0.00005; k2=0.000025; k3=0.0000125
\
\
$R=\frac{||U_k-U_\frac{k}{2}||}{||U_\frac{k}{2}-U_\frac{k}{4}||}=2.025208481097705$

* The order of error is O($h^2$).
檢查圖形是否重疊及R值是否為4
檢驗收斂性 : 分別取N=100,200,400點,time step k=0.01 T=1
\
Then $R=\frac{||U_h-U_\frac{h}{2}||}{||U_\frac{h}{2}-U_\frac{h}{4}||}=4.013069410138102$

## Result

Let $m(t)=\int_{-2}^{2}u(x,t)^2dx$.
The picture is the graph of $\sqrt{m(t)}$

* The numerical solution
N=200; k=0.01 ; total T=20

python code:

import numpy as np
import matplotlib.pyplot as plt

def KSt1(N,k,T):
    # N:the number of points
    # k:time step
    # T:total time
    h=4/N
    r=k/h**4
    s=k/2/h
    alpha=0.05
    
    x=np.arange(-2,2,h)
    
    U=1+0.5*np.exp(-40*x**2)
    
    d0=(6*alpha-2*h**2+(1/r))*np.ones(N)
    D0=np.diag(d0)
    d1=(h**2-4*alpha)*np.ones(N-1)
    D1=np.diag(d1,1)
    D_1=np.diag(d1,-1)
    d2=(alpha)*np.ones(N-2)
    D2=np.diag(d2,2)
    D_2=np.diag(d2,-2)
    D=D0+D1+D_1+D2+D_2
    D[0,N-2]=(alpha)
    D[0,N-1]=(h**2-4*alpha)
    D[1,N-1]=(alpha)

    D[N-2,0]=(alpha)
    D[N-1,0]=(h**2-4*alpha)
    D[N-1,1]=(alpha)
    Ds=r*D
    
    for i in range(round(T/k)):
        U0 = U
        U_old = U0
        
        while np.linalg.norm(Func(U_old,U0,Ds,s))>1e-8:
            g0=np.diag((np.roll(U_old,-1)-np.roll(U_old,1)))# construct the jacobian of nonlinear term
            g_1=-1*np.diag(U_old[1:N],-1)
            g1=1*np.diag(U_old[0:N-1],1)
            g=g0+g1+g_1
            g[0,N-1]=-1*U_old[0]
            g[N-1,0]=1*U_old[N-1]
            g=s*g
            As=Ds-g
            U_new = U_old-np.dot(np.linalg.inv(As),Func(U_old,U0,Ds,s))
            U_old = U_new
        U=U_new
    return U

def Func(Us,U0,D,s):
    F=s*Us*(np.roll(Us,-1)-np.roll(Us,1))
    G=np.dot(D,Us)-F-U0
    return G