---
title: 3. Diffusion equation
tags: Finite difference method
---

# Diffusion equation

**Textbook**
1. Finite Difference Methods for Ordinary and Partial Differential Equations by *Randall J. LeVeque*

---

### Exercise 1

Consider the linear advection-diffusion equation
$$
u_t = u_{xx}-u_x, \quad x\in[0, 100], \quad t\ge 0,
$$
with initial and boundary conditions
$$
u(x,0) = e^{-(x-10)^2}, \quad u(0, t) = u(100,t)=0.
$$
Solve the equation using Crank-Nicolson method.
#### Solution
Given $N\in \mathbf{N}$. Consider the grid points $(x_i,t_n) = (ih,nk),\,i=0,\cdots,N$, where $h=\frac{100}{N+1}$ (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)$.

Using Crank-Nicolson method, we obtain
$$D_t U_i^{n+0.5} = \frac{1}{2} [(D_x^2 U_i^{n+1} -D_x U_i^{n+1}) + (D_x^2 U_i^n -D_x U_i^n)]$$
where
$D_t$ denote the $1^{st}$-order central difference w.r.t. $t$ ,
$D_x^2$ denote the $2^{nd}$-order central difference w.r.t. $x$,
and $D_x$ denote the $1^{st}$-order central difference w.r.t. $x$.

\begin{align*}
\implies \frac{ U_i^{n+1} - U_i^n }{k} &= \frac{1}{2} [\frac{U_{i+1}^{n+1}-2U_i^{n+1}+U_{i-1}^{n+1} }{h^2}-\frac{U_{i+1}^{n+1}-U_{i-1}^{n+1}}{2h}]\\
&+ \frac{1}{2} [\frac{U_{i+1}^n-2U_i^n+U_{i-1}^n }{h^2}-\frac{U_{i+1}^n-U_{i-1}^n}{2h}]\\
\end{align*}
where $i=1,\cdots,N,\, n\geq 0$.

Let $r = \frac{k}{2h^2},\, s = \frac{k}{4h}$. then it can be rewritten as
\begin{align*}
&(s-r) U_{i+1}^{n+1} + (1+2r) U_i^{n+1} - (r+s) U_{i-1}^{n+1}\\
=& (r-s) U_{i+1}^n + (1-2r) U_i^n + (r+s) U_{i-1}^n
\end{align*}
where $i=1,\cdots,N,\, n\geq 0$.

In particular, at $i=1$
$$ (s-r) U_2^{n+1} + (1+2r) U_1^{n+1} = (r-s) U_2^n + (1-2r) U_1^n;$$
at $i=N$
$$ (1+2r) U_N^{n+1} - (r+s) U_{N-1}^{n+1} = (1-2r) U_N^n + (r+s) U_{N-1}^n $$

In matrix form, this is
\begin{align*}
&\begin{pmatrix}
1+2r & s-r & & & \\
-r-s & 1+2r & s-r & & \\
&& \ddots &&\\
&&& \ddots & s-r\\
&&& -r-s & 1+2r
\end{pmatrix}
\begin{pmatrix}
U_1^{n+1} \\ U_2^{n+1} \\ \vdots \\ \vdots \\ U_N^{n+1}
\end{pmatrix}\\
&=\begin{pmatrix}
1-2r & r-s & & & \\
r+s & 1-2r & r-s & & \\
&& \ddots &&\\
&&& \ddots & r-s\\
&&& r+s & 1-2r
\end{pmatrix}
\begin{pmatrix}
U_1^n \\ U_2^n \\ \vdots \\ \vdots \\ U_N^n
\end{pmatrix}
\end{align*}

#### Code(python)

#### Input ####
# grid_size: number of grid points excluding leftmost point x=0.
# 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 LADsolve(grid_size,dt,T): 
    N = grid_size
    k = dt
    # h: mesh width / spacial grid size. dx
    h = 100/N
    
    # iter_max: maximal iteration of time t. iter_max*k < T
    iter_max = int(T/k)
    
    # initial condition u(x,0) = exp(-(x-10)^2)
    x = np.arange(0,100+h,h)
    u = np.zeros( ( N+1 , iter_max+1 ) )
    u[:,0] = np.exp(-(x-10)**2)
    
    r = k/2/h**2
    s = k/4/h
    A = diags([ -r-s , 1+2*r , s-r] , [-1,0,1],shape=(N-1,N-1),format='csc')
    B = diags([ r+s , 1-2*r , r-s ] , [-1,0,1],shape=(N-1,N-1),format='csc')
    
    for i in range(iter_max):
        u[1:N,i+1] = spsolve(A,B*u[1:N,i])
    
    return x,u

#### Result


#### Local Truncation Error
Suppose that $u(x,t)$ is the true solution of our linear advection-diffusion equation, then the local truncation error is
\begin{align*}
&\tau(x_i,t_{n+0.5}) \\
=&
\frac{ u(x_i,t_{n+1}) - u(x_i,t_n) }{k} \\
& - \frac{1}{2}[ \frac{u(x_{i+1},t_{n+1})-2u(x_i,t_{n+1})+u(x_{i-1},t_{n+1}) }{h^2} + \frac{u(x_{i+1},t_n)-2u(x_i,t_n)+u(x_{i-1},t_n) }{h^2} ]\\
& + \frac{1}{2} [\frac{u(x_{i+1},t_{n+1})-u(x_{i-1},t_{n+1})}{2h} + \frac{u(x_{i+1},t_n)-u(x_{i-1},t_n}{2h}]\\
=&
( u_t + \frac{ u_{ttt} }{24} k^3 + \dots ) - ( u_{xx} + \frac{ u_{xxtt} }{8} k^2 + \frac{ u_{xxxx} }{12} h^2 + \dots) + ( u_x + \frac{ u_{xtt} }{8} k^2 + \frac{ u_{xxx} }{6} h^2 + \dots )\\
=&
( \frac{ u_{ttt} }{24} - \frac{ u_{xxtt} }{8} + \frac{ u_{xtt} }{8} ) k^2 + ( -\frac{ u_{xxxx} }{12} + \frac{ u_{xxx} }{6} ) h^2 + o(h^2+k^2),\quad \text{at } (x,t)=(x_i,t_{n+0.5})
\end{align*}

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

1. Consider $N=200$ ( $h = 100/N$) , $k = 0.01$ and $t=100$.

$$R = \frac{\lVert U_h^k - U_{h/2}^k \rVert_2}{\lVert U_{h/2}^k - U_{h/4}^k \rVert_2}=4.060143594084393 \approx 2^2$$

Hence, if time step $k$ is small enough, $\lVert E_h^k\rVert = O(h^2)$, where h is mesh width.
1. Consider $N=10000$ $(h=0.01)$, $k = 0.1$ and $t=100$, then

$$ R = \frac{\lVert U_h^k - U_h^{k/2} \rVert}{\lVert U_h^{k/2} - U_h^{k/4} \rVert} =4.000489890896513 \approx 2^2$$

Hence, if grid size $h$ is small enough, $\lVert E_h^k\rVert = O(k^2)$, where k is time step.
---

### Exercise 2

Consider the linear diffusion equation
$$
u_t = u_{xx}, \quad x\in[-1, 1], \quad t\ge 0,
$$
with initial and boundary conditions
$$
u(x,0) = \sin(\pi x), \quad u(-1, t) = u(1,t)=0.
$$
Solve the equation using Leap-Frog scheme given as
$$
U^{n+2}_i = U^n_i + \frac{2k}{h^2}(U^{n+1}_{i-1}-2U^{n+1}_i+U^{n+1}_{i+1}).
$$

**Solution.**
First devide the interval $[-1,1]$ into $N+1$ pieces, so we insert $N$ points in $(-1,1)$.

* Partition:
$x_0=-1, x_{N+1}=1,$
$x_i=-1+ih$ with $h=\dfrac{2}{N+1}$ , $i=1,...,N$
$t_n = nk$ , $k = \Delta t, n \geq 0$
$U_{i}^{n} = u(x_{i},t_{n})$
$U_{i}^{0} = sin(\pi x_i), U_{0}^{n} = U_{N+1}^{n}=0$
* Unknows: $U_i^{n}=u(x_i,t_n)$ for $i=1,...,N$, $n > 0$

For the first time step, we use Crank-Nicolson method.

Since
$$-rU_{i-1}^1 + (1+2r)U_i^1 - rU_{i+1}^1 = rU_{i-1}^0 + (1-2r)U_i^0 + rU_{i+1}^0 $$
where $$r = \frac{k}{2h^2}, i=1,...,N $$

The corresponding system is
$$
\begin{pmatrix}
(1+r) & -r & & & \\
-r & (1+r) & -r & & \\
& -r & (1+r) & -r & \\
& & & \ddots & \\
& & & -r & (1+r)\\
\end{pmatrix}
\begin{pmatrix}
U_1^1 \\ U_2^1 \\ U_3^1 \\ \vdots \\ \vdots \\ U_N^1
\end{pmatrix}=
\begin{pmatrix}
(1-2r)U_1^0+rU_2^0\\rU_1^0+(1-2r)U_2^0+rU_3^0\\ rU_2^0+(1-2r)U_3^0+rU_4^0\\ \vdots \\rU_{N-1}^0+(1-2r)U_N^0
\end{pmatrix}.
$$

To Solve this system, we use Thomas's algorithm.

Then we use the given Leap-Frog scheme to solve the following time step.

Code:(C++)

#include <iostream>
#include <cmath>
#include <vector>
using namespace std;

const double SMALL = 1.0E-30;

bool tdma( const vector<long double> &a, 
           const vector<long double> &b,
           const vector<long double> &c,
           const vector<long double> &d, 
           vector<long double> &X )
/*
Solving a tri-diagonal system by using the Thomas Algorithm (TDMA).              
a_i X_i-1 + b_i X_i + c_i X_i+1 = d_i,     i = 0, n - 1                                       
Since this is the n x n matrix equation,
a[i], b[i], c[i] are the non-zero diagonals of the matrix and d[i] is the RHS. 
so a[0] and c[n-1] aren't used.                       

Code Refrences:
http://www.cplusplus.com/forum/beginner/247330/                            
*/
{
   int n = d.size();
   vector<long double> P( n, 0 );
   vector<long double> Q( n, 0 );
   X = P;

   // Forward pass
   int i = 0;
   double denominator = b[i];
   P[i] = -c[i] / denominator;
   Q[i] =  d[i] / denominator;
   for ( i = 1; i < n; i++ )
   {
      denominator = b[i] + a[i] * P[i-1];
      if ( abs( denominator ) < SMALL ) return false;
      P[i] =  -c[i]                   / denominator;
      Q[i] = ( d[i] - a[i] * Q[i-1] ) / denominator;
   }

   // Backward pass
   X[n-1] = Q[n-1];
   for (i=n-2;i>=0;i--) X[i]=P[i]*X[i+1]+Q[i];
    return true;
}

int main(void) {
    cout.precision(10);           //More precisely 
    int N;
    double k;                      //N : number of unknown points, k : time step
    cout << "Please enter N(an interger):";
    cin >> N;
    cout << "Please enter k :";
    cin >> k;
    cout << "Please enter t :";
    cin >> t;

    int M = (int) ((t/k)+0.5);
    long double h = 2.0/(N+1);
    long double r = k/(2*pow(h,2.0));
    long double u[M+1][N+2];           //Store u
    vector<long double> x(N+2, 0.0); //Create a vector for x_i
    vector<long double> F(N, 0.0);   //Create a vector for F_i
    vector<long double> U(N, 0.0);   //Create a vector for U_i

    for(int i=0;i<=N+1;i++){         //Devide [-1,1] into N+1 pieces such that x(0) = -1 , x(N+1) = 1
	x[i] = -1 + i*h;             //unknown points are x(1) ~ x(N)
    }
    for(int i=0; i<=N+1; i++)        //initializing u
    {
        u[0][i] = sinl(M_PI*x[i]);
    }

    for(int i=1;i<=N;i++){           //Compute r.h.s of the system
	F[i-1] = r*u[0][i+1] + (1-2*r)*u[0][i] + r*u[0][i-1];
    }

    vector<long double> MainDiag(N,1+2*r);     //Creating the tridiagonal matrix
    vector<long double> LowDiag(N,-r);
    vector<long double> HighDiag(N,-r);
    LowDiag[0] = 0;
    HighDiag[N-1] = 0;
                                          //Compute the first step by Crank-Nicolson method
    tdma(LowDiag,MainDiag,HighDiag,F,U);  //Solve the tridiagonal system by using Thomas's algorithm
	
    for(int i=1; i<=N; i++)              //Update u
    {
        u[1][i] = U[i-1];
    }

    u[1][0] = 0;
    u[1][N+1] = 0;

    int n;
    int n_1;
    int n_2;

    for(int t=2; t<=M; t++)                //Leap-frog
    {
        n = t;
        n_1 = t-1;
        n_2 = t-2;

        for(int i=1; i<=N; i++)              
        {
            u[n][i] = u[n_2][i] + 4*r*(u[n_1][i-1] - 2*u[n_1][i] + u[n_1][i+1]);
        }
        u[n][0] = 0;
        u[n][N+1] = 0;
    }

    cout << "[";               //output u[M+1][N+2]
    for(int i=0; i<=M; i++)
    {
        for(int j=0; j<=N; j++)
        {
            cout << u[i][j] << ",";
        }
        if(i == M) cout << u[i][N+1] << "]";
        else cout << u[i][N+1] << ";";
    }
    
}

Output: matrix of `U[M+1][N+2]`.

The exact solution is $sin(\pi x)exp^{-\pi^2t}$.


**Order:**
Set $N=9,k=0.01$

Set $N=19,k=0.005$

Set $N=39,k=0.0025$

Compute $Err=||u_{exact}-U||_2$ at $t=0.04$ with $N=9,19,39$ and $k=0.01,0.005,0.0025$ respectively.

| $N$ | $h$ | $k$ |$Err$ |
| ---------- |:-------------- |:---------- |:---------- |
| 9 | $0.2000$ | $0.0100$ | $0.009164282576198$ |
| 19 | $0.1000$ | $0.0050$ | $0.002300993600481$ |
| 39 | $0.0500$ | $0.0025$ | $0.000578771135722$ |

Since
$$\frac{Err_{N=19}}{Err_{N=9}}=0.251082785951808 \approx 0.25$$
$$\frac{Err_{N=39}}{Err_{N=19}}=0.251530962798487 \approx 0.25$$
the trucation error is $O(h^2) + O(k^2).$

**Convergence:**
Set $N=9,19,39$ and $k=0.005$ at $t=0.05$

Set $N=19$, $k=0.01,0.005,0.0025$ at $t=0.08$

---

### Exercise 3

Consider the linear diffusion equation
$$
u_t = u_{xx}, \quad x\in[-1, 1], \quad t\ge 0,
$$
with initial and boundary conditions
$$
u(x,0) = \cos(\pi x), \quad u'(-1, t) = u'(1,t)=0.
$$

* Solve the equation using Backward Euler method.
* Show that the mass is conserved in the equation and verify whether or not your scheme preserve the mass.

#### Answer:
* Grid:
To solve the Neumann problem, consider the ghost point method and design the grid
$x_{-1}=-1-h, x_0=-1, ..., x_N=1, x_{N+1}=1+h, h=\frac2{N}$
* Scheme:
The BE method iteratively solves (hence implicitly) $U^{n+1}_:$ using data from $U^n_:$.
(write $r=k/h^2$)
* For interior points, $\frac{U^{n+1}_i-U^n_i}{k}=\frac{U^{n+1}_{i+1}-2U^{n+1}_i+U^{n+1}_{i-1}}{h^2}$
thus $U_i^n = (1+2r)U^{n+1}_i- rU^{n+1}_{i-1}- rU^{n+1}_{i+1}$
* For $i=0$, model the Neumann condition using central difference $\frac{U^n_1-U^n_{-1}}{2h}=0$ and we get the relation $U_0^n = (1+2r)U^{n+1}_0- 2rU^{n+1}_{1}$
* For $i=N$, model the Neumann condition using central difference $\frac{U^{n+1}_{N+1}-U^{n+1}_{N-1}}{2h}=0$ and we get $U^n_N=-2rU^{n+1}_{N-1}+(1+2r)U^{n+1}_N$.

To summarize, we get a `(N+1)x(N+1)` system $$\begin{pmatrix} \frac{1+2r}2 & -r & 0 & ... \\ -r & 1+2r & -r & ... \\ & \ddots & \\ & &1+2r &-r \\ & ...& -r & \frac{1+2r}2 \end{pmatrix}_{N+1\times N+1}\begin{pmatrix} U^{n+1}_0 \\ U^{n+1}_1\\ ... \\ ...\\ U^{n+1}_{N-1} \\ U^{n+1}_{N}\end{pmatrix}_{N+1}=\begin{pmatrix} U^{n}_0/2 \\ U^{n}_1\\ ... \\ ...\\ U^{n}_{N-1} \\ U^{n}_{N}/2\end{pmatrix}_{N+1}$$

* Computation:
We compute `3x3` times with different h's and k's and compare them with the exact solution. $$u(x,t) = e^{-\pi^2t}\cos(\pi x)$$

Python codes -

def headtail(vec):
    n = np.shape(vec)[0]
    veca = np.copy(vec)
    veca[0] = vec[0]/2
    veca[n-1] = vec[n-1]/2
    return veca

def solve_heat_Neumann(n_x, n_t, k):
    grid = np.linspace(-1,1,n_x+1)
    h = grid[1]-grid[0]
    
    F = np.zeros(n_x+1)
    for i in range(n_x+1):
        F[i] = math.cos(math.pi*grid[i])
    
    lamb_t_U = [] #this list collects U(t^n)
    lamb_t_U.append(F)
    
    r = k/h**2
    b = (1+2*r)*np.ones(n_x+1) #diagonal
    b[0] = (1+2*r)/2
    b[n_x] = (1+2*r)/2
    c = (-r)*np.ones(n_x) #superdiagonal
    a = (-r)*np.ones(n_x) #subdiagonal

    
    for i in range(n_t):
        U_new = Thomas(a,b,c, headtail(lamb_t_U[i]))
        lamb_t_U.append(U_new)
    
    return lamb_t_U

Courtesy to cbllelei's implementation of Thomas' algorithm.
https://gist.github.com/cbellei/8ab3ab8551b8dfc8b081c518ccd9ada9


* Error:

Theoretically, $\bigg\{\begin{array}{l}\frac{u(x,t)-u(x,t-k)}k=u_t-1/2 ku_{tt}+O(k^2)\\ \frac{u(x+h,t)-2u(x,t)+u(x-h,t)}{h^2}=u_{xx} +1/6 h^2 u_{xxxx}+O(h^4)\end{array}$
thus $\frac{u(x,t)-u(x,t-k)}{k}-\frac{u(x+h,t)-2u(x,t)+u(x-h,t)}{h^2}=u_t-\frac12ku_{tt} +O(k^2)-u_{xx}-\frac16h^2u_{xxxx}-O(h^4)$

$=-(1/2k - 1/6 h^2)u_{tt}+ O(k^2+h^4)=O(k+h^2)$

Experimentally, I tried 9 combinations and computed their `R`(grid refinement order) and `E`(error from exact solution).

As h refines, R is roughly $0.24$, confirming that this method is $O(h^2)$
; while as k refines, R is roughly $0.3427$, where $-\log_2 0.3427=1.6228$, better than expected $O(k)$.

While fixing $\mu$, (compare diagonally)


`R_1 = 0.3355674678742637`
`R_2 = 0.33774295514543007`
`R_inf = 0.3428220815840617`

which is similar to refining k.

* Mass conservation:

The mass of a system is defined as $\mu(t^*)=\int_{-1}^1 u(x,t^*) dx$. We can use trapezoidal rule to calculate it.

This shows the mass with various h's and k's.
Blue: h=1/64, Orange: h=1/32, Green: h=1/16.


* Codes:

def headtail(vec):
    n = np.shape(vec)[0]
    veca = np.copy(vec)
    veca[0] = vec[0]/2
    veca[n-1] = vec[n-1]/2
    return veca

def solve_heat_Neumann(n_x, n_t, k):
    grid = np.linspace(-1,1,n_x+1)
    h = grid[1]-grid[0]
    
    F = np.zeros(n_x+1)
    for i in range(n_x+1):
        F[i] = math.cos(math.pi*grid[i])
    
    lamb_t_U = [] #this list collects U(t^n)
    lamb_t_U.append(F)
    
    r = k/h**2
    b = (1+2*r)*np.ones(n_x+1) #diagonal
    b[0] = (1+2*r)/2
    b[n_x] = (1+2*r)/2
    c = (-r)*np.ones(n_x) #superdiagonal
    #c[0] = -2*r
    a = (-r)*np.ones(n_x) #subdiagonal
    #a[n_x-1] = -2*r
    
    realmat = np.diag(a,-1)+np.diag(b,0)+np.diag(c,1)
    
    for i in range(n_t):
        U_new = Thomas(a,b,c, headtail(lamb_t_U[i]))
        lamb_t_U.append(U_new)
    
    return lamb_t_U

Courtesy to cbllelei's implementation of Thomas' algorithm.
https://gist.github.com/cbellei/8ab3ab8551b8dfc8b081c518ccd9ada9

---

### Exercise 4

Consider the diffusion equation
$$
u_t = u_{xx}-e^{\sin(x)}, \quad x\in[0, 1], \quad t\ge 0,
$$
with initial and boundary conditions
$$
u(x,0) = \cos(\pi x), \quad u'(0, t) =1, \quad u(1,t)=0.
$$
Solve the equation and find the long-time behavior of the solution.
#### Answer:
* Partition:
$x_0=0, \ x_{N+1}=1,\ x_i=ih$ with $h=\dfrac{1}{N+1},\ i=1,...,N$ \
k is the time step , $t_n=nk$
$U_{N+1}^n=u(x_{N+1},t_n)=u(1,t_n)=0.$

* Unknown: $U_i^n=u(x_i,t_n)$ for $i=0,1,...,N$

$u_t = u_{xx}-e^{\sin(x)},
\quad\Rightarrow\quad\dfrac{U_{i}^{n+1}-U_{i}^{n}}{k}= \dfrac{U_{i+1}^{n+1}-2U_{i}^{n+1}+U_{i-1}^{n+1}}{h^2}-e^{\sin(x_i)}$

For $i=1,2,\dots ,N-1$
let $r=\frac{k}{h^2}$ and we have $$(1+2r)U_i^{n+1}-rU_{i+1}^{n+1}-rU_{i-1}^{n+1}=U_i^{n}-ke^{\sin(x_i)}$$

For i=N
$$(1+2r)U_N^{n+1}-rU_{N+1}^{n+1}-rU_{N-1}^{n+1}=U_N^{n}-ke^{\sin(x_N)}$$
Since $U_{N+1}^{n+1}=0$ ,we get $(1+2r)U_N^{n+1}-rU_{N-1}^{n+1}=U_N^{n}-ke^{\sin(x_N)}$

* For i=0 find a ghost point $U_{-1}^{n+1}$ so that we can approximate $u'(0,t)$

Since $u'(0,t)=1\approx\dfrac{U_1^{n+1}-U_{-1}^{n+1}}{2h}\quad\Rightarrow\quad U_{-1}^{n+1}=U_1^{n+1}-2h$.

For i=0
$$(1+2r)U_0^{n+1}-rU_{1}^{n+1}-rU_{-1}^{n+1}=U_0^{n}-ke^{\sin(x_0)}$$

Substituting $U_{-1}^{n+1}=U_1^{n+1}-2h$.

we get $(1+2r)U_0^{n+1}-2rU_{1}^{n+1}=U_0^{n}-ke^{\sin(x_0)}-2rh$

Hence the system becomes
$$
\begin{pmatrix}
\frac{1+2r}{2}& -r & & \\
-r & 1+2r & -r & & \\
& -r &1+2r & -r & \\
& &\ddots & \ddots \\
& & & -r & 1+2r\\
\end{pmatrix}
\begin{pmatrix}
U_0^{n+1} \\ U_1^{n+1} \\ \vdots \\ \vdots \\ U_N^{n+1}
\end{pmatrix}=
\begin{pmatrix}
\frac{U_0^{n}-ke^{\sin(x_0)}-2rh}{2}\\U_1^{n}-ke^{\sin(x_1)}\\U_2^{n}-ke^{\sin(x_2)}\\ \vdots \\U_N^{n}-ke^{\sin(x_N)}
\end{pmatrix}.
$$

$u_t-u_{xx}+e^{\sin(x_i)}-\dfrac{U_{0}^{n+1}-U_{0}^{n}}{k}+ \dfrac{U_{1}^{n+1}-2U_{0}^{n+1}-U_{i-1}}{h^2}-e^{\sin(x_i)}=O(h^2)+O(k)$

檢驗收斂性 : 分別取N= 99,199,399點,time step k=0.01 疊代次數 n=500
$h\approx\frac{1}{N+1}$ Then $\frac{||U_\frac{h}{2}-U_\frac{h}{4}||}{||U_h-U_\frac{h}{2}||}=0.24999847521995003\approx\frac{1}{4}$

Python code:

def heat(N,k,n):
    #N : the number of points
    #k : time step
    #n : iterative times
    import numpy as np

    h=1/(N+1)
    r=k/(h**2)
    x=np.linspace(0,1,N+2)
    U0=np.cos(np.pi*x[0:N+1]) #U0 initial
    U=U0+(-1)*k*np.exp(np.sin(x[0:N+1]))
    U[0]=0.5*(U[0]-2*h*r)

    d0=(1+2*r)*np.ones(N+1)
    d0[0]=0.5*(1+2*r)
    D0=np.diag(d0)

    d1=-r*np.ones(N)
    D1=np.diag(d1,1)

    d2=-r*np.ones(N)
    D2=np.diag(d2,-1) 
    D=D0+D1+D2

    Dinv=np.linalg.inv(D)

    for i in range(n):
        U=np.dot(Dinv,U)
        U=U+(-1)*k*np.exp(np.sin(x[0:N+1]))
        U[0]=0.5*(U[0]-2*h*r)

    import matplotlib.pyplot as plt
    #plt.plot(range(1,U.size),U[1:U.size])
    plt.plot(x[1:N+1],U[1:U.size])
    return U

the long-time behavior of the solution.

---

### Exercise 5

Consider the nonlinear diffusion equation
$$
u_t = u_{xx}, \quad x\in[0, 1], \quad t\ge 0,
$$
with initial and boundary conditions
$$
u(x,0) = \cos(\pi x/2), \quad u'(0, t) =0, \quad u(1,t)=\sin(t).
$$
Solve the equation and find the long-time behavior of the solution.

**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 = u_{xx}, \quad x\in[0, 1], \quad t\ge 0, \quad
u(x,0) = \cos(\pi x/2), \quad u'(0, t) =0, \quad u'(1,t)=0.
$$
exact solution:
$$u(x,t)=C_0+\sum_{k=1}^{\infty} C_k e^{-k^2 \pi^2 t} \cos (k\pi x) \quad , C_0=\frac{2}{\pi} \quad , C_k=2\int_{0}^{1} u(x,0) \cos (k \pi x) dx$$
We use central difference to approximate $u_{xx}$
\
$$
\frac{1}{h^2}(U_{i+1}-2U_i+U_{i-1})=u_{xx}
$$
With partition(ghost point):
$\mu=0.16=\frac{h^2}{k^2}$ (fixed)
$k=0.1, 0.05, 0.025$ (time step)
$N=\frac{1}{h}$ and $h=\sqrt{0.16} \ k=0.04, 0.02, 0.01$
$x_i=ih-0.5h$ with $i=0 \sim N+1$
$x_{0}\approx x_1, \\ x_{N+1}\approx x_N,$
(With truncation error $O(h^2)$)

then we get the system : $\frac{dU}{dt}= A\,U$
$$\frac{dU}{dt}=\frac{1}{h^2} \begin{pmatrix}
-1 & 1 & 0 & \cdots & 0 \\
1 & -2 & 1 & \cdots & 0 \\
\vdots & \ddots & \ddots & \ddots & \vdots \\
0 & \cdots & 1 & -2 & 1 \\
0 & \cdots & 0 & 1 & -1
\end{pmatrix}U
\\
\\
\text{with } U(x_i,0)=U^0_i=
\begin{pmatrix}
\cos (\frac{\pi x_1}{2})\\
\vdots\\
\cos (\frac{\pi x_i}{2})\\
\vdots\\
\cos (\frac{\pi x_N}{2})
\end{pmatrix}
$$

Now we get an ode problem,
so we can use differential equation solvers of julia to solve it.
\
Method: Midpoint - The second order midpoint method
(With truncation error $O(k^2)$)

mu=0.16 ##=h^2/k^2
k1=0.1
h1=mu^0.5*k1
N1=Int64(round(1/h1, digits=0))
x1=(1:N1).*h1.-(0.5*h1);
T=4*pi;
du_c = repeat([1],N1-1);
d_c = [-1;repeat([-2],N1-2);-1];
dl_c = repeat([1],N1-1);
A_c = Tridiagonal{Float64}(dl_c,d_c,du_c)
u0 = cos.(pi.*x1/2);
f(u,p,t) = A_c*u/h1^2
tspan = (0.0,T)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob, Midpoint(), reltol=1e-8, abstol=1e-8,save_everystep=false,dt=k1)


because we know the exact solution of this problem is
$$u(x,t)=C_0+\sum_{k=1}^{\infty} C_k e^{-k^2 \pi^2 t} \cos k\pi x \quad , C_0=\frac{2}{\pi} \quad , C_k=2\int_{0}^{1} u(x,0) \cos (k \pi x) dx$$
we can get the error of long time behavior$(t=4 \pi)$

|$h=\sqrt{\mu}k$ | k=0.1 | k=0.05 | k=0.025 |
| -------- | -------- | -------- |-------- |
|1 norm|0.000104732 | 2.61807e-5 | 6.54503e-6|
|2 norm|0.000104732 | 2.61807e-5 | 6.54503e-6|
|$\infty$ norm|0.000104736 | 2.61837e-5 | 6.54825e-6|
which is $O(k^2)=O(h^2)$

#### True case:
Now let $u(1,t)=\sin(t)$
approximate $U_{N+1}$:
(ghost point with truncation error $O(h^2)$)
$$
\sin(t) = u(1,t) \approx \frac{U_{N}+U_{N+1}}{2} \Rightarrow U_{N+1} \approx 2 \sin(t)- U_N \\
$$
then we get the system : $\frac{dU}{dt}= A\,U+p(t)$

$$\frac{dU}{dt}=\frac{1}{h^2} \begin{pmatrix}
-1 & 1 & 0 & \cdots & 0 \\
1 & -2 & 1 & \cdots & 0 \\
\vdots & \ddots & \ddots & \ddots & \vdots \\
0 & \cdots & 1 & -2 & 1 \\
0 & \cdots & 0 & 1 & -3
\end{pmatrix}U+\frac{1}{h^2}
\begin{pmatrix}
0\\
0\\
\vdots\\
0\\
2 \sin t
\end{pmatrix}
\\
\\
\text{with } U(x_i,0)=U^0_i=
\begin{pmatrix}
\cos (\frac{\pi x_1}{2})\\
\vdots\\
\cos (\frac{\pi x_i}{2})\\
\vdots\\
\cos (\frac{\pi x_N}{2})
\end{pmatrix}
$$


#### long time behavior of the solution:

which looks like a periodic solution

---