notebook.community

Edit and run



In [ ]:

    
#https://nbviewer.jupyter.org/github/waltherg/notebooks/blob/master/2013-12-03-Crank_Nicolson.ipynb



In [10]:

    
import numpy
from matplotlib import pyplot









    



/Users/saigerutherford/anaconda2/lib/python2.7/site-packages/matplotlib/font_manager.py:280: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.
  'Matplotlib is building the font cache using fc-list. '



In [11]:

    
numpy.set_printoptions(precision=3)



In [12]:

    
L = 1.
J = 100
dx = float(L)/float(J-1)
x_grid = numpy.array([j*dx for j in range(J)])



In [13]:

    
T = 200
N = 1000
dt = float(T)/float(N-1)
t_grid = numpy.array([n*dt for n in range(N)])



In [14]:

    
D_v = float(10.)/float(100.)
D_u = 0.01 * D_v

k0 = 0.067
f = lambda u, v: dt*(v*(k0 + float(u*u)/float(1. + u*u)) - u)
g = lambda u, v: -f(u,v)
 
sigma_u = float(D_u*dt)/float((2.*dx*dx))
sigma_v = float(D_v*dt)/float((2.*dx*dx))

total_protein = 2.26



In [16]:

    
no_high = 10
U =  numpy.array([0.1 for i in range(no_high,J)] + [2. for i in range(0,no_high)])
V = numpy.array([float(total_protein-dx*sum(U))/float(J*dx) for i in range(0,J)])



In [20]:

    
ylim=((0., 2.1))
xlabel=('x'); ylabel=('concentration')
pyplot.plot(x_grid, U)
pyplot.plot(x_grid, V)
pyplot.show()



In [21]:

    
A_u = numpy.diagflat([-sigma_u for i in range(J-1)], -1) +\
      numpy.diagflat([1.+sigma_u]+[1.+2.*sigma_u for i in range(J-2)]+[1.+sigma_u]) +\
      numpy.diagflat([-sigma_u for i in range(J-1)], 1)
        
B_u = numpy.diagflat([sigma_u for i in range(J-1)], -1) +\
      numpy.diagflat([1.-sigma_u]+[1.-2.*sigma_u for i in range(J-2)]+[1.-sigma_u]) +\
      numpy.diagflat([sigma_u for i in range(J-1)], 1)
        
A_v = numpy.diagflat([-sigma_v for i in range(J-1)], -1) +\
      numpy.diagflat([1.+sigma_v]+[1.+2.*sigma_v for i in range(J-2)]+[1.+sigma_v]) +\
      numpy.diagflat([-sigma_v for i in range(J-1)], 1)
        
B_v = numpy.diagflat([sigma_v for i in range(J-1)], -1) +\
      numpy.diagflat([1.-sigma_v]+[1.-2.*sigma_v for i in range(J-2)]+[1.-sigma_v]) +\
      numpy.diagflat([sigma_v for i in range(J-1)], 1)



In [22]:

    
print A_u









    



[[ 1.981 -0.981  0.    ...,  0.     0.     0.   ]
 [-0.981  2.962 -0.981 ...,  0.     0.     0.   ]
 [ 0.    -0.981  2.962 ...,  0.     0.     0.   ]
 ..., 
 [ 0.     0.     0.    ...,  2.962 -0.981  0.   ]
 [ 0.     0.     0.    ..., -0.981  2.962 -0.981]
 [ 0.     0.     0.    ...,  0.    -0.981  1.981]]



In [23]:

    
f_vec = lambda U, V: numpy.multiply(dt, numpy.subtract(numpy.multiply(V, 
                     numpy.add(k0, numpy.divide(numpy.multiply(U,U), numpy.add(1., numpy.multiply(U,U))))), U))



In [24]:

    
print f(U[0], V[0])









    



0.00996135898275



In [25]:

    
print f(U[-1], V[-1])









    



-0.0623832232232



In [26]:

    
print f_vec(U, V)









    



[ 0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01
  0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01
  0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01
  0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01
  0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01
  0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01
  0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01
  0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01   0.01
  0.01   0.01  -0.062 -0.062 -0.062 -0.062 -0.062 -0.062 -0.062 -0.062
 -0.062 -0.062]



In [27]:

    
U_record = []
V_record = []

U_record.append(U)
V_record.append(V)

for ti in range(1,N):
    U_new = numpy.linalg.solve(A_u, B_u.dot(U) + f_vec(U,V))
    V_new = numpy.linalg.solve(A_v, B_v.dot(V) - f_vec(U,V))
    
    U = U_new
    V = V_new
    
    U_record.append(U)
    V_record.append(V)



In [28]:

    
ylim=((0., 2.1))
xlabel=('x'); ylabel=('concentration')
pyplot.plot(x_grid, U)
pyplot.plot(x_grid, V)
pyplot.show()



In [36]:

    
U_record = numpy.array(U_record)
V_record = numpy.array(V_record)

fig, ax= subplots()
xlabel=('x'); ylabel=('t')
heatmap = ax.pcolor(x_grid, t_grid, U_record, vmin=0., vmax=1.2)









    



---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-36-486b16ddffc2> in <module>()
      2 V_record = numpy.array(V_record)
      3 
----> 4 fig, ax= subplots()
      5 xlabel=('x'); ylabel=('t')
      6 heatmap = ax.pcolor(x_grid, t_grid, U_record, vmin=0., vmax=1.2)

NameError: name 'subplots' is not defined



In [ ]: