The diffusion equations:

Assuming a constant diffusion coefficient, D, we use the Crank-Nicolson methos (second order accurate in time and space):

u[n+1,j]-u[n,j] = 0.5 a {(u[n+1,j+1] – 2u[n+1,j] + u[n+1,j-1])+(u[n,j+1] – 2u[n,j] + u[n,j-1])}

A linear system of equations, A.U[n+1] = B.U[n], should be solved in each time setp.

As an example, we take a Gaussian pulse and study variation of density with time. Such example can occur in several fields of physics, e.g., quantum mechanics.

As stated above, we have two arrays: A and B. We form them once, and then calculate inverse(A).B:

#———————————————–

# populate the coefficient arrays

#———————————————-

from scipy.linalg import svd

from scipy.linalg import inv

for j in arange(0, nx_step):

coefa[j, (j-1)%nx_step] = – lamb

coefa[j, (j)%nx_step] = 2*lamb + 2.

coefa[j, (j+1)%nx_step] = – lamb

#———————————————–

lamb = – lamb

for j in arange(0, nx_step):

coefb[j, (j-1)%nx_step] = – lamb

coefb[j, (j)%nx_step] = 2*lamb + 2.

coefb[j, (j+1)%nx_step] = – lamb

#———————————————–

coefa = scipy.linalg.inv(coefa) # inverse of A

coef = numpy.dot(coefa, coefb)

coef = scipy.linalg.inv(coef)

for i in arange(1,nt_step): #———– main loop ———

ans = scipy.linalg.solve(coef, uo)

uo = ans

plot(uoz,’k’, ans,’r’)

draw()

t[i] = t[i-1] + dt

print dt*nt_step – t[i], ‘ ‘, ans.sum()

rho[:,i] = ans

The result is shown in the following figure: as expected, the peak of the density profile falls down. If one continues the simulation for a longer time, the density distribution will be completely uniform.