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=D∂2ρ(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,j−2ρi,j+ρi−1,jΔx2⇒ρi,j+1=(1−2λ)ρi,j+λ(ρi+1,j+ρi−1,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:
*Sovling an example system with initial value of ρcenter=1.0 at center and boundary conditions with ρL=ρR=0.
A3. Implicit scheme (backward time)
In implicit scheme, time evolves backwards.
ρi,j−ρi,j−1Δt=Dρi+1,j−2ρi,j+ρi−1,jΔx2⇒(1+2λ)ρi,j−λ(ρi+1,j+ρi−1,j)=ρi,j−1
where λ=D⋅ΔtΔx2.
Rewrite this in matrix form,
A→ρi,j−R=→ρi,j−1
Supplemented by initial condition: ρi,0=ρ(xi),
and Boundary conditions: ρ0,j=ρ0, ρm,j=ρm, we can solve it like this:
*Sovling an example system with initial value of ρcenter=1.0 at center and boundary conditions with ρL=ρR=0.
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,j−2ρi,j+ρi−1,jΔx2+ρi+1,j+1−2ρi,j+1+ρi−1,j+1Δx2)⇒(2+2λ)ρi,j+1−λ(ρi+1,j+1+ρi−1,j+1)=(2−2λ)ρi,j+λ(ρi+1,j+ρi−1,j)
where λ=D⋅ΔtΔx2
Rewrite this in matrix form,
A→ρi,j+1−R=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:
Solvlng an example system with initial value of ρcenter=1.0 at center and boundary conditions with ρL=ρR=0.
B. Solving diffusion equation with a potential profile.
∂ρ(x,t)∂t=D∂∂x[∂ρ(x,t)∂x+1kBT∂U(x)∂xρ(x,t)]=D∂2ρ(x,t)∂x2+DkBT⋅U′(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,j−2ρi,j+ρi−1,jΔx2+ρi+1,j+1−2ρi,j+1+ρi−1,j+1Δx2)+D⋅U′(xi)2kBT(ρi+1,j+1−ρi−1,j+12Δx+ρi+1,j−ρi−1,j2Δx)+D⋅U′′(xi)2kBT(ρi,j+1+ρi,j)⇒(2+2λ+τi)ρi,j+1+(ηi−λ)ρi+1,j+1−(λ+ηi)ρi−1,j+1=(2−2λ−τi)ρi,j+(λ−ηi)ρi+1,j+(λ+ηi)ρi−1,j
where λ=D⋅ΔtΔx2, ηi=−D⋅Δt2kBTΔx⋅U′(xi), τi=−D⋅ΔtkBT⋅U′′(xi).
Rewrite this in matrix form,
A→ρi,j+1−R=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:
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 used in this example and its derivatives are as follows:
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=D∂∂x[∂ρ0(x,t)∂x+1kBT∂U0(x)∂xρ0(x,t)]−[k0→1+kwf0(x)]ρ0(x,t)+k1→0ρ1(x,t)∂ρ1(x,t)∂t=D∂∂x[∂ρ1(x,t)∂x+1kBT∂U1(x)∂xρ1(x,t)]−[k1→2+k1→0+kwf1(x)]ρ1(x,t)+k2→1ρ2(x,t)+k0→1ρ0(x,t)∂ρ2(x,t)∂t=D∂∂x[∂ρ2(x,t)∂x+1kBT∂U2(x)∂xρ2(x,t)]−[k2→3+k2→1+kwf2(x)]ρ2(x,t)+k3→2ρ3(x,t)+k1→2ρ1(x,t)⋮⋮∂ρn−1(x,t)∂t=D∂∂x[∂ρn−1(x,t)∂x+1kBT∂Un−1(x)∂xρn−1(x,t)]−[kn−1→n+kn−1→n−2+kwfn−1(x)]ρn−1(x,t)+kn→n−1ρn(x,t)+kn−2→n−1ρn−2(x,t)∂ρn(x,t)∂t=D∂∂x[∂ρn(x,t)∂x+1kBT∂Un(x)∂xρn(x,t)]−[kn→n−1+kwfn(x)]ρn(x,t)+kn−1→nρn−1(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,j−2ρki,j+ρi−1,jΔx2+ρi+1,j+1−2ρi,j+1+ρi−1,j+1Δx2)+D⋅U′k(xi)2kBT(ρki+1,j+1−ρki−1,j+12Δx+ρki+1,j−ρki−1,j2Δx)+D⋅U′′k(xi)2kBT(ρki,j+1+ρki,j)−kk→k−1+kk→k+1+kwfk(xi)2(ρki,j+1+ρki,j)+kk−1→k2(ρk−1i,j+1+ρk−1i,j)+kk+1→k2(ρk+1i,j+1+ρk+1i,j)⇒(2+2λ+τki)ρki,j+1+(ηi−λ)ρki+1,j+1−(λ+ηki)ρki−1,j+1+ξkiρki,j+1−ζkρk+1i,j+1−κkρk−1i,j+1=(2−2λ−τki)ρki,j+(λ−ηki)ρki+1,j+(λ+ηki)ρki−1,j−ξkiρki,j+1+ζkρk+1i,j+1+κkρk−1i,j+1
where
λ=D⋅ΔtΔx2ηki=−D⋅Δt2kBTΔx⋅U′k(xi)τki=−D⋅ΔtkBT⋅U′′k(xi)ξki=(kk→k+1+kk→k−1+kwfk(xi))Δtζk=kk+1→kΔtκk=kk−1→kΔt
We further define following (m−1)×(m−1) matrices:
Rewrite in matrix form,
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.
The above example calculations are produced using Julia 6.0 with matplotlib. You can find the example jupyter notebook file HERE.