Notes on numerically solving diffusion equations.

2018/1/7 14:46 PM posted in  Numerical 0 Comments

Many phenomena in various science fields are mathematically expressed by using the well-known evolution equations. The diffusion equation is one of them, and mathematically corresponds to the Markov process in relation to the normal distribution rule.

A. Standard diffusion equation

Standard diffusion equation is a simplistic second order partial differential equation written as:

ρ(x,t)t=D2ρ(x,t)x2

A1. Finite difference method for time and space

Finite difference of time:

ρt(xi,tj)=ρ(xi,tj+Δt)ρ(xi,tj)Δt+O(Δt)

Finite difference of space:

2ρx2(xi,tj)=ρ(xi+Δx,tj)2ρ(xi,tj)+ρ(xiΔx,tj)dx2+O(Δx2)

A2. Explicit scheme (forward time)

Apply equation(2) and (3) to equation(1), and let ρi,j=ρ(xi,tj), we get the expression of forward time evolution

 ρi,j+1ρi,jΔt=Dρi+1,j2ρi,j+ρi1,jΔx2ρi,j+1=(12λ)ρi,j+λ(ρi+1,j+ρi1,j)

where λ=DΔtΔx2

.
Rewrite this in matrix form,

ρi,j+1=Aρi,j+R(accounts for boundary conditions)

Supplemented by initial condition: ρi,0=ρ(xi),
and Boundary conditions: ρ0,j=ρ0, ρm,j=ρm, we can solve it like this:

Screenshot 2018-01-07 15.04.40

*Sovling an example system with initial value of ρcenter=1.0 at center and boundary conditions with ρL=ρR=0.

forward

A3. Implicit scheme (backward time)

In implicit scheme, time evolves backwards.

 ρi,jρi,j1Δt=Dρi+1,j2ρi,j+ρi1,jΔx2(1+2λ)ρi,jλ(ρi+1,j+ρi1,j)=ρi,j1

where λ=DΔtΔx2.


Rewrite this in matrix form,

Aρi,jR=ρi,j1

Supplemented by initial condition: ρi,0=ρ(xi),
and Boundary conditions: ρ0,j=ρ0, ρm,j=ρm, we can solve it like this:

Screenshot 2018-01-07 15.18.52

*Sovling an example system with initial value of ρcenter=1.0 at center and boundary conditions with ρL=ρR=0.

backward

A4. Crank-Nicolson scheme

In Crank-Nicolson scheme, we take the average of explicit scheme and implicit scheme.

 ρi,j+1ρi,jΔt=D2(ρi+1,j2ρi,j+ρi1,jΔx2+ρi+1,j+12ρi,j+1+ρi1,j+1Δx2)(2+2λ)ρi,j+1λ(ρi+1,j+1+ρi1,j+1)=(22λ)ρi,j+λ(ρi+1,j+ρi1,j)

where λ=DΔtΔx2

.
Rewrite this in matrix form,

Aρi,j+1R=Bρi,j+R

Supplemented by initial condition: ρi,0=ρ(xi),
and Boundary conditions: ρ0,j=ρ0, ρm,j=ρm, we can solve it like this:

Screenshot 2018-01-07 15.25.12

Solvlng an example system with initial value of ρcenter=1.0 at center and boundary conditions with ρL=ρR=0.

crank-nicolson

B. Solving diffusion equation with a potential profile.

ρ(x,t)t=Dx[ρ(x,t)x+1kBTU(x)xρ(x,t)]=D2ρ(x,t)x2+DkBTU(x)ρ(x,t)x+DkBTU(x)ρ(x,t)

The finite difference in Crank-Nicolson scheme for the above equation is:

 ρi,j+1ρi,jΔt=D2(ρi+1,j2ρi,j+ρi1,jΔx2+ρi+1,j+12ρi,j+1+ρi1,j+1Δx2)+DU(xi)2kBT(ρi+1,j+1ρi1,j+12Δx+ρi+1,jρi1,j2Δx)+DU(xi)2kBT(ρi,j+1+ρi,j)(2+2λ+τi)ρi,j+1+(ηiλ)ρi+1,j+1(λ+ηi)ρi1,j+1=(22λτi)ρi,j+(ληi)ρi+1,j+(λ+ηi)ρi1,j

where λ=DΔtΔx2, ηi=DΔt2kBTΔxU(xi), τi=DΔtkBTU(xi).

Rewrite this in matrix form,

Aρi,j+1R=Bρi,j+R

Supplemented by initial condition: ρi,0=ρ(xi),
and Boundary conditions: ρ0,j=ρ0, ρm,j=ρm, we can solve it like this:

Screenshot 2018-01-07 15.28.00

Solving an example system with initial value of ρ=1.0 at x=30 and boundary conditions with ρL=ρR=0 using the potential profile of nw=0.

potential-crank-nicolson

Potential used in this example and its derivatives are as follows:

Screenshot 2018-01-07 15.33.08
Screenshot 2018-01-07 15.33.23
Screenshot 2018-01-07 15.33.34

C. Solving coupled diffusion equations with potential profiles and reactions.

When reaction terms are added to the equation(16), equations corresponding to different nw become coupled equations:

{ρ0(x,t)t=Dx[ρ0(x,t)x+1kBTU0(x)xρ0(x,t)][k01+kwf0(x)]ρ0(x,t)+k10ρ1(x,t)ρ1(x,t)t=Dx[ρ1(x,t)x+1kBTU1(x)xρ1(x,t)][k12+k10+kwf1(x)]ρ1(x,t)+k21ρ2(x,t)+k01ρ0(x,t)ρ2(x,t)t=Dx[ρ2(x,t)x+1kBTU2(x)xρ2(x,t)][k23+k21+kwf2(x)]ρ2(x,t)+k32ρ3(x,t)+k12ρ1(x,t)ρn1(x,t)t=Dx[ρn1(x,t)x+1kBTUn1(x)xρn1(x,t)][kn1n+kn1n2+kwfn1(x)]ρn1(x,t)+knn1ρn(x,t)+kn2n1ρn2(x,t)ρn(x,t)t=Dx[ρn(x,t)x+1kBTUn(x)xρn(x,t)][knn1+kwfn(x)]ρn(x,t)+kn1nρn1(x,t)

Let ρki,j=ρk(xi,tj), the finite difference in Crank-Nicolson scheme for nw=k is:

 ρki,j+1ρki,jΔt=D2(ρki+1,j2ρki,j+ρi1,jΔx2+ρi+1,j+12ρi,j+1+ρi1,j+1Δx2)+DUk(xi)2kBT(ρki+1,j+1ρki1,j+12Δx+ρki+1,jρki1,j2Δx)+DUk(xi)2kBT(ρki,j+1+ρki,j)kkk1+kkk+1+kwfk(xi)2(ρki,j+1+ρki,j)+kk1k2(ρk1i,j+1+ρk1i,j)+kk+1k2(ρk+1i,j+1+ρk+1i,j)(2+2λ+τki)ρki,j+1+(ηiλ)ρki+1,j+1(λ+ηki)ρki1,j+1+ξkiρki,j+1ζkρk+1i,j+1κkρk1i,j+1=(22λτki)ρki,j+(ληki)ρki+1,j+(λ+ηki)ρki1,jξkiρki,j+1+ζkρk+1i,j+1+κkρk1i,j+1

where

λ=DΔtΔx2ηki=DΔt2kBTΔxUk(xi)τki=DΔtkBTUk(xi)ξki=(kkk+1+kkk1+kwfk(xi))Δtζk=kk+1kΔtκk=kk1kΔt

We further define following (m1)×(m1) matrices:

Screenshot 2018-01-07 15.38.20

Rewrite in matrix form,

Screenshot 2018-01-07 15.38.29

Solving an example system with initial value of ρ=1.0 at x=40 for all nw=0:9 and boundary conditions with ρL=ρR=0 using the potential profile of nw=0.
coupled

The above example calculations are produced using Julia 6.0 with matplotlib. You can find the example jupyter notebook file HERE.