---
title: 5. Fisher-KPP equation
tags: Finite difference method
---

---

### Problem 3 (亞衛)

@RyanHsieh 1. 你的方法以及你的 code 在收斂階數上有沒有一致需要 check. 2. 所以原題在 $t=40$ 的時候的解究竟長怎樣? 在你的回答裡我看到三條不一樣的線, 所以是誰? 還是都不是?

@RyanHsieh 1. stability 一開始假設的有點跟我不一樣, 沒有$\lambda$? 2. Von Neumann 應該只能做線性的部分, 非線性那一項應該做不動. 3. 看起來 x=0 的邊界有錯, 解在那邊怪怪的.

Fisher-KPP equation
$$
u_t = u_{xx}+u-u^2, \quad x\in[0, 100], \quad t\in[0, 40],
$$
with initial and boundary conditions
$$
u(x,0) = e^{-x}, \quad u(0, t) =1, \quad u(100,t)=0.
$$

#### **solution**.
We fisrt consider the same equation with boundary condition:
$u(x,0) = (1+\exp(\sqrt{\dfrac{1}{6}}x))^{-2}, \\
u(0,t) = (1+\exp(\dfrac{-5}{6}t))^{-2}, \text{and} \\
u(200,t) = (1+\exp(\dfrac{-5}{6}t+\sqrt{\dfrac{1}{6}}\times 200))^{-2}.$
The Exact solution is
$$u(x,t)=(1+\exp(\dfrac{-5}{6}t+\sqrt{\dfrac{1}{6}}x))^{-2}.
$$
First devide the interval $[0,100]$ into $N+1$ pieces, so we insert $N$ points in $(0,100)$.
Nonation : Here we use $U_i^n$ instead of using $u(x_i,t_n)$ ,the numerical approximation.

* **Partition**:
$x_0=0, \quad x_N=100,\quad x_i=ih$ with $h=\dfrac{100}{N+1}$ for $i=0,1,2,...,N,N+1$.
$t_0=0, \quad t_M=40,\quad t_n=nk$ with $k=\dfrac{40}{M}$ for $n=0,1,...,M$.
* **Unknows**: $U_i^n=u(x_i,t_n)$ for $i=0,1,...,N-1$ and $n=1,2,...,M$.

#### Central Difference Method
Here we use the Backward Euler method, so the original equation
$$u_t = u_{xx}+u-u^2$$ becomes
$\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}
+U_i^n(1-U_i^n)$.
Hence
$-rU_{i+1}^{n+1}+(1+2r)U_i^{n+1}-rU_{i-1}^{n+1}
= (1+k-kU_i^n)U_i^n$ where $r=\dfrac{k}{h^2}$.

So the correponding system is
$$
\begin{pmatrix}
1+2r & -r & & & \\
-r & 1+2r & -r & & \\
& -r & 1+2r & -r & \\
& & & \ddots & \\
& & & -r & 1+2r\\
\end{pmatrix}
\begin{pmatrix}
U_1^{n+1} \\ U_2^{n+1} \\ \vdots \\ \vdots \\ U_{N-1}^{n+1}
\end{pmatrix}=
\begin{pmatrix}
(1+k-kU_1^n)U_1^{n}+rU^{n+1}_0 \\ (1+k-kU_2^n)U_2^{n} \\ \vdots \\ \vdots \\ (1+k-kU_{N-1}^n)U_{N-1}^{n}+rU^{n+1}_N
\end{pmatrix}\quad\quad(*)
$$

#### Truncation error
$\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}
+U_i^n(1-U_i^n)$

The local truncation error is
\begin{align*}
\tau(x_i,t_n) & =D_tu_i^n-D_{xx}u_i^{n+1}-u_i^{n}(1-u_i^n) \\
& =(ku_t+\frac{k^2}{2}u_{tt}+\frac{k^3}{6}u_{ttt}+\dots)-(u_{xx}+\frac{h^2}{12}u_{xxxx}+\dots)-u(1-u)\\
&=(u_t-u_{xx}-u(1-u)) + (k-1)u_t+\frac{k^2}{2}u_{tt}+\dots-\frac{h^2}{12}u_{xxxx}+\dots \\
&\approx O(h^2) + O(k).
\end{align*}

---

__2020.07.02 Update__

##### Matlab code

clc;clear;

% Initailization
N=10;
h=100/N;
k=0.1;
r=k/(h^2);
T=40;
M=T/k;
xx=[0:h:100];
tt=[0:k:T];

% Boundary Condition
Sol=zeros(M+1,N+1);
v=zeros(N-1,1);
Sol(1,:) = (1+exp(sqrt(1/6)*xx)).^(-2);
Sol(:,1) = (1+exp((-5/6)*tt)).^(-2);
Sol(:,N+1) = (1+exp((-5/6)*tt+sqrt(1/6)*100)).^(-2);

% Tridiagonal system
S = diag((1+2*r)*ones(1,N-1))+diag((-r)*ones(1,N-2),1)+diag((-r)*ones(1,N-2),-1);
x_new = Sol(1,:);

% Central Difference Method
for i=1:M
    x_old = x_new;
    v(1) = r*Sol(i+1,1);
    v(N-1) = r*Sol(i+1,N+1);
    right = (((1+k-k.*Sol(i,2:N)).*Sol(i,2:N))' + v);
    x_new = S\right;
    Sol(i+1,:)=[Sol(i+1,1),(x_new)',Sol(i+1,N+1)];
end

% Creating Exact solution
ts = (-5/6)*tt';
xs = sqrt(1/6)*xx;
ET = repmat(ts,1,N+1);
EX = repmat(xs,M+1,1);
Exact = (1+exp(ET+EX)).^(-2);

% Calculate the Error
Error = abs(Exact - Sol);
err = norm(Error,inf)

Setting $T=40$
| $N$ | $h$ | $k$ |$k/h^2$|$Err$ | $\frac{Err_h}{Err_{2h}}$|
| :---------- |:---|:-------------- |:---------- |:---------- |:-|
| $10$ | $10$ | $0.1$ | $10^{-3}$ | $3.195845119677438$| |
| $20$ | $5$ | $0.025$ | $10^{-3}$ | $3.754365127925062$| $1.174764416713660$|
| $40$ | $2.5$ | $0.00625$ | $10^{-3}$ | $2.869063707583049$|$0.764194107345309$ |
| $80$ | $1.25$ | $0.0015625$ | $10^{-3}$ | $1.593838088400404$|$0.555525513144873$ |
| $160$ | $0.625$ | $0.000390625$ | $10^{-3}$ | $0.815406831116279$| $0.511599538905882$|
| $320$ | $0.3125$ | $0.00009765625$ | $10^{-3}$ | $0.409762213268130$|$0.502524871795803$ |
Exact Solution when $T=40$ with $N=80(h=1.25)$ and $k=0.0015625$:
Numerical Solution when $T=40$ with $N=80(h=1.25)$ and $k=0.0015625$:

Setting $T=1$
| $N$ | $h$ | $k$ |$k/h^2$|$Err$ | $\frac{Err_h}{Err_{2h}}$|
| :---------- |:---|:-------------- |:---------- |:---------- |:-|
| $10$ | $10$ | $0.1$ | $10^{-3}$ | $0.004976076044017$| |
| $20$ | $5$ | $0.1$ | $10^{-3}$ | $0.002158588593650$|$0.433793329232857$ |
| $40$ | $2.5$ | $0.1$ | $10^{-3}$ | $0.005017586842186$| $2.324475750935783$|
| $80$ | $1.25$ | $0.1$ | $10^{-3}$ | $0.002908354883233$| $0.579632196652908$|
| $160$ | $0.625$ | $0.1$ | $10^{-3}$ | $0.001483586242644$| $0.510111833737019$|
| $320$ | $0.3125$ | $0.1$ | $10^{-3}$ | $0.000746273196146$|$0.503019760291128$ |
| $640$ | $0.078125$ | $0.1$ | $10^{-3}$ | $0.000373472669036$|$0.500450332351122$ |

Exact Solution when $T=1$ with $N=80(h=1.25)$ and $k=0.0015625$:
Numerical Solution when $T=1$ with $N=80(h=1.25)$ and $k=0.0015625$:

Convergence:
$T=1$
$T=10$
$T=40$

From the table above, we can check that the trucation error is $O(h^2) + O(k).$

To solve the original equation, the code is the same except

Sol(1,:) = exp(-xx);
Sol(:,1) = 1;
Sol(:,N+1) = 0;

Result:
$T=40$ with $N=160(h=0.625)$ and $k=0.15625$.
Convergence:
$N=10,20,40,80,160$ at $T=40$

References:
1. [Numerical Study of One Dimensional Fishers KPP Equation with Finite Difference Schemes.](https://pdfs.semanticscholar.org/1263/aed1db64f0ea7577e6d22f15930aa379a064.pdf)
2. [A Numerical Treatment of Fisher Equation.](https://cyberleninka.org/article/n/611770.pdf)

---

### Problem 4 (沂君)

@7G35H_3DTJyS0go4-C366Q good job, 不過 0.26 -> 0.25 -> 0.22 看起來不像會收斂到 0.25 的樣子? N 再往上取試試. 如果要跑很久, 算到 T=1 就好.

Fisher-KPP equation
$$
u_t = u_{xx}+u-u^2, \quad x\in[0, 200], \quad t\in[0, 40],
$$
with initial and boundary conditions
$$
u(x,0) = e^{-x}, \quad u(0, t) =1, \quad u(200,t)=0.
$$

#### **Solution.**
First we solve the equation with different initial and boundary condition to check our code.

$$
u_t = u_{xx}+u-u^2, \quad x\in[0, 200], \quad t\in[0, 40],
$$

$$
u(x,0) = (1+exp(\sqrt{\frac{1}{6}}x))^{-2}, \quad u(0, t) =(1+exp(-\frac{5}{6}t))^{-2},
$$

$$
u(200,t)=(1+exp(-\frac{5}{6}t+\sqrt{\frac{1}{6}}\times 200))^{-2}.
$$

Exact solution:
$$
u(x,t)=(1+exp(-\frac{5}{6}t+\sqrt{\frac{1}{6}} x))^{-2}
$$

Devide the interval $[0,200]$ into $N+1$ pieces, so we insert $N$ points in $(0,200)$.

* Partition:
$x_0=0, x_{N+1}=200,$
$x_i=0+ih$ with $h=\dfrac{200}{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} = (1+exp(\sqrt{\frac{1}{6}}x_i))^{-2}$
* Initial and boundary conditions:
$U_{0}^{t} = (1+exp(-\frac{5}{6}t))^{-2}$
$U_{200}^{t} = (1+exp(-\frac{5}{6}t+\sqrt{\frac{1}{6}}\times 200))^{-2}$
* Unknows: $U_i^{n}=u(x_i,t_n)$ for $i=1,...,N$, $n > 0$

Using Crank-Nicolson method, we obtain
$$
\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-2U_i^n+U_{i-1}^n }{h^2}]\\
+ \frac{U_i^{n+1}+U_i^n}{2} - (\frac{U_i^{n+1}+U_i^n}{2})^2
$$

$$
\implies (-\frac{1}{2}r)U_{i-1}^{n+1}+(1+r-\frac{k}{2})U_i^{n+1}+(-\frac{1}{2}r)U_{i+1}^{n+1}\\
+(-\frac{1}{2}r)U_{i-1}^{n}+(-1+r-\frac{k}{2})U_i^n+(-\frac{1}{2}r)U_{i+1}^n\\
+k(\frac{U_i^{n+1}+U_i^n}{2})^2=0=F(U), \quad r =\frac{k}{h^2}
$$

In each time step $t=(n+1)k$, we use newton's method to find $U_i^{n+1}$.

That is, find $U^m= \left[\begin{matrix} U_1^{n+1} \\ U_2^{n+1} \\ \vdots \\ U_N^{n+1} \\\end{matrix}\right]$ by $U^{m+1} = U^m - (J[U^m])^{-1}F[U^m]$

with initial guass $U^0=0$, $J$ is Jocabian matrix of $F$ and consider $U_i^n$ as known constant.

So we get the system

$U^{m+1}=$
$$
U^m-
\begin{pmatrix}
1+r-\frac{k}{2}+\frac{k}{2}(U_1^{n+1}+U_1^{n}) & -\frac{1}{2}r & & & \\
-\frac{1}{2}r & 1+r-\frac{k}{2}+\frac{k}{2}(U_2^{n+1}+U_2^{n}) & -\frac{1}{2}r & & \\
& & & \ddots & \\
\end{pmatrix}^{-1}
F(U^m).
$$

##### Pseudo code
1. $X = \{0+ih | i=0,\dots,N+1,\, h=\frac{200}{N+1} \},\quad M=\frac{40}{k}$
1. $U_i^0 = (1+exp(\sqrt{\frac{1}{6}} x_i))^{-2}$
1. $\text{for }t=1,\dots,M\\
\quad U^0=[(1+exp(-\frac{5}{6}t))^{-2} \quad 0\dots0 \quad(1+exp(-\frac{5}{6}t+\sqrt{\frac{1}{6}} \times 200))^{-2}]^t\\
\quad \text{while not convergent:}\\
\quad \quad U^k \leftarrow U^{k-1}\\
\quad \quad U^k \leftarrow U^k - J(U^k)^{-1} F(U^k)\\
\quad \text{endwhile}\\
\quad U_i^t=U^k(i-1), \quad i=1,\dots,N\\
\quad U_0^t=(1+exp(-\frac{5}{6}t))^{-2}\\
\quad U_{200}^t=(1+exp(-\frac{5}{6}t+\sqrt{\frac{1}{6}} \times 200))^{-2}\\
\text{endfor}$

##### Code:(C++)

#include <iostream>
#include <cmath>
#include <vector>
#include <fstream>
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;
}

long double infinity_norm(const vector<long double> &a) {
    long double max=0;
    long double k;
    for(int i = 0; i < a.size(); i++)
    {
        k = fabs(a[i]);
        if(k > max) max = k;
    } 
    return max;
}

int main(void) {

    ofstream outfile("./u.txt");

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

    int M = (int) ((40/k)+0.5);       //max t = 40
    long double h = 200.0/(N+1);
    long double r = k/(pow(h,2.0));
    long double u[2][N+2]={0};           //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> V(N, 0.0);   //Create a vector for V_i

    for(int i=0;i<=N+1;i++){         //Devide [0,200] into N+1 pieces such that x_0 = 0 , x_(N+1) = 200
    	x[i] = 0 + i*h;             //unknown points are x_1 ~ x_(N)
    }
    
    for(int i=0; i<=N+1; i++)        //initializing u
    {
        u[0][i] = pow( ( 1 + ( exp( pow(0.166666666666667,0.5)*x[i]))) ,-2); 
    }
    
    if(outfile.is_open()){
        outfile << "[";
        for(int j=0; j<=N; j++)
        {
            outfile << u[0][j] << ",";
        }
        if(0 == M) outfile << u[0][N+1] << "]";
        else outfile << u[0][N+1] << ";";   
    }
    

    int n;
    int n_1;

    for(int t=1; t<=M; t++)                
    {
        n = t%2;                            //Time n
        n_1 = ((n%2)+2-1)%2;                //Time n-1
        
        u[n][0] = pow( 1+exp( (-0.833333333333333)*t*k ) ,-2);
        u[n][N+1] = pow( 1+exp( (-0.833333333333333)*t*k + pow(0.166666666666667,0.5)*200) ,-2);
        
        long double Err = 1;
        long double Tol = 1.0E-16;
        vector<long double> Unew(N,0.0);
        vector<long double> Uold_Unew(N);

        while(Err > Tol)
        {            
            for(int i=0; i<N; i++)      //update u
            {
                u[n][i+1] = Unew[i];
            }

            for(int i=1;i<=N;i++){     
	        F[i-1] = (-0.5)*r*u[n][i+1] + (1+r-0.5*k)*u[n][i] + (-0.5)*r*u[n][i-1] 
                        + (-0.5)*r*u[n_1][i+1] + (-1+r-0.5*k)*u[n_1][i] + (-0.5)*r*u[n_1][i-1]
                        + k*pow( (u[n][i]+u[n_1][i])/2, 2) ;
	    }

            vector<long double> MainDiag(N,1+r-0.5*k);     //Creating the tridiagonal matrix J
	    vector<long double> LowDiag(N,-0.5*r);
	    vector<long double> HighDiag(N,-0.5*r);
            LowDiag[0] = 0;
            HighDiag[N-1] = 0;
            for(int i=0; i<N; i++)
            {
                MainDiag[i] = MainDiag[i] + 0.5*k*u[n][i] + 0.5*k*u[n_1][i];
            }
        
            tdma(LowDiag,MainDiag,HighDiag,F,V);

            //V = J^(-1)F
            for(int i=0; i<N; i++)
            {
                Unew[i] = u[n][i+1]-V[i];
            }

            for(int i=0; i<N; i++)
            {
                Uold_Unew[i] = u[n][i+1]-Unew[i];
            }

            Err = infinity_norm(Uold_Unew);
        }

        for(int i=0; i<N; i++)              //Update u
        {
            u[n][i+1] = Unew[i];
        }
        
        if(outfile.is_open()){
            //if(t==M){ outfile << "[";
            for(int j=0; j<=N; j++)
            {
                outfile << u[n][j] << ",";
            }
            if(t == M) outfile << u[n][N+1] << "]";
            else outfile << u[n][N+1] << ";";
            //}
        }
        
    }
}

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

Exact solution:

Set $N=199, k=0.1$

At t=40.

Compute $Err=||u_{exact}-U||_1$ at $t=40$ with $N=99,199,399,799,1599,3199,6399$ and $k=0.2,0.1,0.05,0.025,0.0125,0.00625,0.003125$ respectively.

| $N$ | $h$ | $k$ |$Err$ | $\frac{Err_h}{Err_{2h}}$|
| ---------- |:-------------- |:---------- |:---------- |:-|
| 99 | $2.00000$ | $0.200000$ | $5.563056214616273$ | |
| 199 | $1.00000$ | $0.100000$ | $1.466021177852882$ |$0.263528017926744$ |
| 399 | $0.50000$ | $0.050000$ | $0.370390248021738$ |$0.252649998251872$ |
| 799 | $0.25000$ | $0.025000$ | $0.092811838248455$ |$0.250578514807462$ |
| 1599 | $0.12500$ | $0.012500$ | $0.023217880269299$ |$0.250160762974497$ |
| 3199 | $0.06250$ | $0.006250$ | $0.005808800758313$ |$0.250186523960759$ |
| 6399 | $0.03125$ | $0.003125$ | $0.001456418463975$ |$0.250726186793426$ |

From the table above, we can check that the trucation error is $O(h^2) + O(k^2).$

To solve the original equation, the code is the same except

for(int i=0; i<=N+1; i++)        //initializing u  //line 89
{
    u[0][i] = exp(-x[i]);
}
u[0][0]=1;
u[0][N+1]=0;

for(int t=1; t<=M; t++)  //line 101              
{
        n = t%2;                            //Time n
        n_1 = ((n%2)+2-1)%2;                //Time n-1
        
        u[n][0] = 1;         
        u[n][N+1] = 0;

Set $N=199, k=0.1$

Convergence:

---