0
$\begingroup$

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?

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

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.subplots_adjust(left=0.18,bottom=0.15)
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X,T,sum_mat,cmap=cm.jet,lw=0.1,rstride=1, cstride=1)
plt.xlabel("$x$")
plt.ylabel("$t$")
ax.set_zlabel('$u$')
fig.savefig('ans.png',transparent=True,dpi=300)

Finally, this link is the plot of the numerical solution.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.