# Solving the diffusion equation with state-dependent diffusion coefficient numerically?

I am having an issue numerically solving the following diffusion equation with state-dependent diffusion coefficient and need some help finding out what to do. The equation I wish to solve is:

$$\frac{\partial u}{\partial t} = \frac{\partial}{\partial x} \left( D(x) \frac{\partial u}{\partial x} \right)$$

where $D(x)$ is some sigmoidal function.

I am struggling to solve the equation analytically but am trying to use the Crank-Nicholson finite difference method to implement it numerically. Can anybody provide some insight into analytical/numerical techniques available to solve this equation?

• Do you mean $$\frac{\partial u}{\partial t} = \frac{\partial}{\partial x} \left( D(x) \frac{\partial u}{\partial x} \right) ?$$ – Professor Vector Aug 16 '17 at 11:04
• Apologies, yes! – K.Pentland Aug 16 '17 at 11:39
• Can it be assumed that $x \leq 1$? – ekkilop Aug 16 '17 at 16:13
• Yes x can be less than or equal to one if necessary! – K.Pentland Aug 17 '17 at 14:04

This can be solved numerically using implicit Euler method. The partial differential equation:

$$\frac{\partial u}{\partial t}=D(x)\frac{\partial^2 u}{\partial x^2}+\frac{d D(x)}{d x}\frac{\partial u}{\partial x}$$ Subjected to the following initial and boundary conditions: \begin{align*} u(x,0)&=\sin(\pi x)\\ u(0,t)&=0\\ u(1,t)&=0 \end{align*} With finite difference, we can write: $$\frac{u_{i,j+1}-u_{i,j}}{\Delta t}=D(x_i)\big[\frac{u_{i+1,j+1}-2u_{i,j+1}+u_{i-1,j+1}}{\Delta x^2}\big]+\frac{d D(x)}{d x}\Big|_{x=x_i}\cdot\big[\frac{u_{i+1,j+1}-u_{i,j+1}}{\Delta x}\big]$$ Note that $i$ refers to the $i$th node of the discretized space, and $j$ corresponds to the $j$th node of the discretized time. The equation can be rearranged, such that: $$u_{i,j+1}-D(x_i)\big[\frac{u_{i+1,j+1}-2u_{i,j+1}+u_{i-1,j+1}}{\Delta x^2}\big]\Delta t-\frac{d D(x)}{d x}\Big|_{x=x_i}\cdot\big[\frac{u_{i+1,j+1}-u_{i,j+1}}{\Delta x}\big]\Delta t=u_{i,j}$$ Now let us write the terms in the following form: \begin{align*} A_{i,0}&=-D(x_i)\frac{\Delta t}{\Delta x^2}\\ A_{i,1}&=1+2D(x_i)\frac{\Delta t}{\Delta x^2}+\frac{d D(x)}{d x}\Big|_{x=x_i}\frac{\Delta t}{\Delta x}\\ A_{i,2}&=-D(x_i)\frac{\Delta t}{\Delta x^2}-\frac{d D(x)}{d x}\Big|_{x=x_i}\frac{\Delta t}{\Delta x} \end{align*} Therefore, at each time step we solve this equation: $$\begin{bmatrix} A_{1,1} & A_{1,2} & 0 & 0&\cdots & 0 & 0 & 0\\ A_{2,0} & A_{2,1} & A_{2,2} & 0&\cdots & 0 & 0 & 0\\ 0 & A_{3,0} & A_{3,1} & A_{3,2}&\cdots & 0 & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \ddots & \ddots & \vdots& \vdots \\ 0 & 0 & 0 & 0&\cdots& A_{N-3,0} & A_{N-3,1} & A_{N-3,2} \\ 0 & 0 & 0 & 0&\cdots & 0 & A_{N-2,0} & A_{N-2,1} \end{bmatrix}\begin{bmatrix} u_{1,j+1}\\u_{2,j+1}\\u_{3,j+1}\\ \vdots \\u_{N-3,j+1}\\u_{N-2,j+1} \end{bmatrix}=\begin{bmatrix} u_{1,j}\\u_{2,j}\\u_{3,j}\\ \vdots \\u_{N-3,j}\\u_{N-2,j} \end{bmatrix}$$ With $u_0$ and $u_{N-1}$ always equal to zero. Assume $D(x)$ is a sigmoidal function: $$D(x)=\frac{1}{1+e^{-x}}$$ $$\frac{dD(x)}{dx}=\frac{e^{-x}}{(1+e^{-x})^2}$$ We then implement this using Python 3.5. The program is as follows:

import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.cm as cm
def f1(x):
return 1/(1+np.exp(-x))
def f2(x):
return np.exp(-x)/(1+np.exp(-x))**2
dx=0.01
x=np.linspace(0.01,0.99,99)
x2=np.linspace(0,1,101)
dt=0.001
t_step=500
sum_mat=np.zeros((int(t_step/25+1),len(x)+2))
u=np.zeros(len(x))
for i in range(len(x)):
u[i]=np.sin(np.pi*x[i])
t=[]
for j in range(t_step):
if j==0:
t.append(j*dt)
c_count=0
for i in range(len(x)):
sum_mat[c_count][i+1]=u[i]
c_count=c_count+1
else:
Am=np.zeros((len(x),len(x)))
for i in range(len(x)):
if i==0:
Am[i][i]=1+f1(x[i])*2*dt/dx**2+f2(x[i])*dt/dx
Am[i][i+1]=-f1(x[i])*dt/dx**2-f2(x[i])*dt/dx
elif i==len(x)-1:
Am[i][i]=1+f1(x[i])*2*dt/dx**2+f2(x[i])*dt/dx
Am[i][i-1]=-f1(x[i])*dt/dx**2
else:
Am[i][i]=1+f1(x[i])*2*dt/dx**2+f2(x[i])*dt/dx
Am[i][i-1]=-f1(x[i])*dt/dx**2
Am[i][i+1]=-f1(x[i])*dt/dx**2-f2(x[i])*dt/dx
u_next=np.linalg.solve(Am,u)
u=u_next
if (j+1)%25==0:
t.append(j*dt)
for i in range(len(x)):
sum_mat[c_count][i+1]=u_next[i]
u[i]=u_next[i]
c_count=c_count+1
X,T=np.meshgrid(x2,t)
fig = plt.figure()
plt.rc('font', **{'family': 'serif', 'serif': ['Computer Modern']})
plt.rc('text', usetex=True)
plt.gcf().set_size_inches(4,3,forward=True)
plt.xlabel("$x$")
plt.ylabel("$t$")
ax.set_zlabel('$u$')