notebook.community

Edit and run

Heat Equation

\begin{eqnarray*} \frac{\partial u(x,y,t)}{\partial t} & = & a\nabla^2 u(x,y,t) =a\left( \ \frac{\partial^2 u(x,y,t)}{\partial x^2}+\frac{\partial^2 u(x,y,t)}{\partial y^2}\right) \end{eqnarray*}

Numeric Scheme

  • discretization: $u(x_i,y_j,t_m)= U^m_{i,j}$
  • Equations \begin{eqnarray} \frac{\partial u}{\partial t} &\sim& \frac{U^{m+1}{i,j}-U^m{i,j}}{\triangle t}\ \frac{\partial^2 u}{\partial x^2} &\sim& \frac{ U^m{i+1,j}-2U^m{i,j}+U^m{i-1,j}}{\triangle x^2}\ \frac{\partial^2 u}{\partial y^2} &\sim& \frac{ U^m{i,j+1}-2U^m{i,j}+U^m{i,j-1}}{\triangle y^2}\ \nabla^2 u(x,y) &\sim& \frac{ U^m{i+1,j}-2U^m{i,j}+U^m_{i-1,j}}{\triangle x^2}+ \ 
             \frac{ U^m_{i,j+1}-2U^m_{i,j}+U^m_{i,j-1}}{\triangle y^2}
    \end{eqnarray}
  • The simulated data at time, $m+1$, step:
$$ U^{m+1}_{i,j}=U^m_{i,j}+ a\triangle t\left(\ \frac{ U^m_{i+1,j}-2U^m_{i,j}+U^m_{i-1,j}}{\triangle x^2}+\frac{ U^m_{i,j+1}-2U^m_{i,j}+U^m_{i,j-1}}{\triangle y^2}\ \right)$$

Stability

$$\triangle t \le \frac{1}{2a}\frac{(\triangle x \triangle y)^2}{\triangle x^2+\triangle y^2}$$

Example

Initial case: $(x,y)\in[0,1]^2$ $$u(x,y)=\left\{\begin{array}{ll} 1, &0.1\le\text{ for } (x-0.5)^2+(y-0.5)^2\le0.2, \\ 0, &\text{ otherwise } \end{array}\right.$$





In [1]:



    


%pylab inline
import numpy as np
import matplotlib.pyplot as plt
from numpy import exp,sin,pi,tanh
from matplotlib import animation
from JSAnimation import IPython_display
import time



    


















    






Populating the interactive namespace from numpy and matplotlib

















In [2]:



    


dx=0.01        
dy=0.01        
a=0.5          
timesteps=500  

nx = int(1/dx)
ny = int(1/dy)

dx2=dx**2 
dy2=dy**2 
dt = dx2*dy2/( 2*a*(dx2+dy2) )



    











In [48]:



    


ui = np.zeros([nx,ny])
u = np.zeros([nx,ny])

for i in np.arange(nx):
   for j in np.arange(ny):
      if ( ( (i*dx-0.5)**2+(j*dy-0.5)**2 <= 0.2)
         & ((i*dx-0.5)**2+(j*dy-0.5)**2>=.1) ):
            ui[i,j] = 1

def evolve_ts(u, ui):
   u[1:-1, 1:-1] = ui[1:-1, 1:-1] + \
    a*dt*( (ui[2:, 1:-1] - 2*ui[1:-1, 1:-1] + ui[:-2, 1:-1])/dx2 + (ui[1:-1, 2:] - 2*ui[1:-1, 1:-1] + ui[1:-1, :-2])/dy2 )



    











In [4]:



    


timesteps=1
tstart = time.time()

for m in np.arange(1, timesteps+1):
   evolve_ts(u, ui)
   #print "Computing u for m =", m
tfinish = time.time()
print "Done."
print "Total time: ", tfinish-tstart, "s"
print "Average time per time-step using numpy: ", ( tfinish - tstart )/timesteps, "s."



    


















    






Done.
Total time:  0.000778913497925 s
Average time per time-step using numpy:  0.000778913497925 s.

















In [3]:



    


def evolve_ts(u, ui, a,dt,dx,dy,dx2,dy2):
    u[1:-1, 1:-1] = ui[1:-1, 1:-1] + \
       a*dt*( (ui[2:, 1:-1] - 2*ui[1:-1, 1:-1] + ui[:-2, 1:-1])/dx2 + (ui[1:-1, 2:] - 2*ui[1:-1, 1:-1] + ui[1:-1, :-2])/dy2 ) 
def UI(nx,ny,a):
    dx = 1/float(nx)
    dy = 1/float(ny)
    ui = np.zeros([nx,ny])
    for i in np.arange(nx):
        for j in np.arange(ny):
            if ( ( (i*dx-0.5)**2+(j*dy-0.5)**2 <= 0.2) & ((i*dx-0.5)**2+(j*dy-0.5)**2>=.1) ):
               ui[i,j] = 1
    return ui            
def heat_sim(nx,ny,a,UI,timesteps):         
    dx = 1/float(nx)
    dy = 1/float(ny)

    dx2=dx**2 
    dy2=dy**2 
    dt = dx2*dy2/( 2*a*(dx2+dy2) )
    
    u = np.zeros([nx,ny])
    ui=UI
    for m in np.arange(1, timesteps+1):
        evolve_ts(u, ui, a,dt,dx,dy,dx2,dy2)
        ui = np.copy(u)  
    return ui



    











In [52]:



    


nx,ny=100,100
a=0.5
#u = np.zeros([nx,ny])
ui0=UI(nx,ny,a)
ui=heat_sim(nx,ny,a,ui0,1000)
im=plt.imshow(ui,cmap=cm.hot, interpolation='nearest', origin='lower',extent=[0,1,0,1])
im.set_clim(vmin=0, vmax=1)
plt.colorbar()



    


















    
Out[52]:








<matplotlib.colorbar.Colorbar instance at 0x10ca13200>










    
























In [4]:



    


anim = plt.figure(figsize=(6,6))
ax = anim.add_subplot(111)
ax.set_title("2D Heat Equation",fontsize=14)

plt.xlim([0,1])
plt.ylim([0,1])

def init():
    nx,ny=100,100
    a=0.5
    #u = np.zeros([nx,ny])
    ui0=UI(nx,ny,a)
    ui=heat_sim(nx,ny,a,ui0,1)
    return ax.imshow(ui,cmap=cm.hot, interpolation='nearest', origin='lower',extent=[0,1,0,1],vmax=1)

def animate(i):
    nx,ny=100,100
    a=0.5
    #u = np.zeros([nx,ny])
    ui0=UI(nx,ny,a)
    heat_sim(nx,ny,a,ui0,i*50)
    ui=heat_sim(nx,ny,a,ui0,i*100)
    return ax.imshow(ui,cmap=cm.hot, interpolation='nearest', origin='lower',extent=[0,1,0,1],vmax=1)

 
animation.FuncAnimation(anim, animate, init_func=init, frames=20, interval=50)



    


















    
Out[4]:











    
    

    
    

    
    
    
    
    
    
    
    
    
  

     Once 
     Loop 
     Reflect 
  























In [ ]:

Fisher Model

\begin{eqnarray*} \frac{\partial N(x,t)}{\partial t} & = & D\frac{\partial^2 N(x,t)}{\partial x^2} +rN(x,t)\left(1-\frac{N(x,t)}{K}\right) \cr \frac{\partial N }{\partial n}|_{\partial \Omega}&=&0 \\ N(x,0) &=& f(x) \end{eqnarray*}
  • $N=0,1$: equilibrium points.
  • $0<N(x,0)<1$:
    • $c\ge c_{\min}=2$: there are infinit (increasing) travelling wave solution, $V(z)$, \begin{eqnarray*} N(x,t) &=& V(\pm x+ c t)\\ \lim_{z\to-\infty}V(z)&=&1 \\ \lim_{z\to\infty}V(z)&=&=0 \end{eqnarray*}
    • $c=5/\sqrt6$, for constant $C$, $$V(z)=\left[1+C\exp\left(-\frac{z}{\sqrt6}\right)\right]^{-2}$$
    • $0<c<2$: No travelling wave solution

Example

  • $D=1,r=1,K=1$
  • $V(z)=0.5\tanh(z)+0.5$
  • discretization: $N(t_k,x_i)= U^k_i$

  • Equations \begin{eqnarray} \frac{\partial N}{\partial t} &\sim& \frac{U^{k+1}_i-U^k_i}{\triangle t}\\ \frac{\partial^2 N}{\partial x^2} &\sim& \frac{ U^k_{i+1}-2U^k_i+U^k_{i-1}}{\triangle x^2}\\ \end{eqnarray}

  • The simulated data at time, $m+1$, step:

$$ U^{k+1}_{i}=U^k_i+ \triangle t\left(\ r U^k_i(1-U^k_i)+D\frac{ U^k_{i+1}-2U^k_{i}+U^k_{i-1}}{\triangle x^2}\ \right)$$
  • Boundary Condition: \begin{eqnarray} \frac{\partial N }{\partial n}|{\partial \Omega}=0 &\to& N^k{-1}-N^k_0=0, N^kn-N^k{n+1}=0 \ 
      \text{(visual points)}   &\to& N^k_{-1}=N^k_0,N^k_{n+1}= N^k_n \\
    N^{k+1}_0 &=& U^k_0+ \triangle t\left(\ r U^k_0(1-U^k0)+D\frac{ U^k{1}-U^k_{0}}{\triangle x^2} \right)\ N^{k+1}_n &=& U^k_n+ \triangle t\left(\ r U^k_n(1-U^kn)+D\frac{ U^k{n-1}-U^k_{n}}{\triangle x^2} \right) \end{eqnarray}
  • convergent condition: $$ \Delta t\le \frac{\Delta x^2}{4r\Delta x^2+2D}$$




In [2]:



    


z=np.linspace(-10,10,100)
V=lambda z: (1+np.exp(z/2.))**(-2.)
plt.plot(z,V(z))



    


















    
Out[2]:








[<matplotlib.lines.Line2D at 0x105ef5790>]










    
























In [51]:



    


D=0.2
r=1.
K=1.
c=3.

dx=0.1
dt=0.001
b0=-20.;b1=20.;
nx=int((b1-b0)/dx)+1
timesteps=500  
dx2=dx**2



    











In [52]:



    


dx2/(r*dx2+D)



    


















    
Out[52]:








0.04761904761904762

















In [47]:



    


def evolve_ts(u, ui):
   u[1:-1] = ui[1:-1] + \
     dt*( (ui[2:] - 2*ui[1:-1] + ui[:-2])/dx2 +r*ui[1:-1]*(1.-ui[1:-1]))
   u[0]=ui[0]+dt*(D*(ui[1] - ui[0])/dx2 + r*ui[0]*(1.-ui[0]))
   u[-1]=ui[-1]+dt*(D*(ui[-2] - ui[-1])/dx2 + r*ui[-1]*(1.-ui[-1]))     
   ui = np.copy(u)     
   return ui



    











In [87]:



    


ui = np.zeros([nx])
u = np.zeros([nx])
x=np.linspace(-20,20,nx)
f=lambda x,t: (1+np.exp((x+c*t)/2.))**(-2.)

for i in np.arange(nx):
    ui[i]=f(x[i],0)
plt.xlim(-20,20)
plt.ylim(-0.2,1.2)
for i in np.arange(5000):
    ui=evolve_ts(u,ui)
    if (i%1200==0):
       plt.plot(x,ui,'r--')
plt.text(4,1,'Travelling Wave',size=16) 
plt.text(0,0.5, r'$\circ\to$',size=36,color='green')



    


















    
Out[87]:








<matplotlib.text.Text at 0x1116ebe90>










    
























In [5]:



    


D=0.1
r=1.
K=1.
c=1.

def fish_sim(dx,dt,b0,b1,c,UI,timesteps): 
    nx=int((b1-b0)/dx)+1
    dx2=dx**2. 

    ui = UI
    U= np.zeros([nx,timesteps])
    U[:,0]=ui
    
    u = np.zeros([nx])

    def evolve_ts(u, ui):                       
        u[1:-1] = ui[1:-1] + dt*( D*(ui[2:] - 2*ui[1:-1] + ui[:-2])/dx2 + r*ui[1:-1]*(1.-ui[1:-1]) )
        # Boundary Points
        u[0]=ui[0]+dt*(D*(ui[1] - ui[0])/dx2 + r*ui[0]*(1.-ui[0]))
        u[-1]=ui[-1]+dt*(D*(ui[-2] - ui[-1])/dx2 + r*ui[-1]*(1.-ui[-1]))

        ui = np.copy(u)
        return ui

    for m in np.arange(1, timesteps):         
        ui=evolve_ts(u, ui)
        U[:,m]=ui
    return U



    











In [6]:



    


dx=0.1
dt=0.02
b0=-20.;b1=20.;
nx=int((b1-b0)/dx)+1
x=np.linspace(-20,20,nx)
ui = np.zeros([nx])
u = np.zeros([nx])
#dt = 0.00001
f=lambda x,t: 0.5*np.tanh(-(x+c*t))+0.5
f=lambda x,t: (1+np.exp((x+c*t)/2.))**(-2.)

for i in np.arange(nx):
    ui[i]=f(x[i],0)
timesteps=1200



    











In [137]:



    


dx2=dx*dx
dx2/(r*dx2+D)
1/(4*r+2*D/dx2)



    


















    
Out[137]:








0.04166666666666667

















In [127]:



    


dx,nx,timesteps,len(ui)



    


















    
Out[127]:








(0.1, 401, 1200, 401)

















In [7]:



    


F=fish_sim(dx,dt,-20,20,c,ui,timesteps+1)



    











In [136]:



    


plt.xlim(-20,20)
plt.ylim(-.2,1.2)
ttime=int(timesteps*dt)
plt.text(0,0.5,'time %d' %ttime)
plt.plot(x,F[:,-1])



    


















    
Out[136]:








[<matplotlib.lines.Line2D at 0x11db7c610>]










    
























In [8]:



    


anim = plt.figure(figsize=(8,6))

ax = plt.axes(xlim=(-20, 20), ylim=(-0.2, 1.2))

ax.set_title("Fisher Model",fontsize=14)
line,=ax.plot([],[])
frames=100

time_text1 = ax.text(9, 1.1, '',size=14)
time_text2 = ax.text(14, 1.1, '',size=14)


def init():
    time_text1.set_text('time = ')
    line.set_data([],[])
    return line,time_text1,time_text2
    
def animate(i):
    frame=timesteps/frames
    ttext=i*float(frame*dt)
    time_text2.set_text( ttext )
    line.set_data(x,F[:,(i)*frame])
    return line,time_text1,time_text2
 
animation.FuncAnimation(anim, animate, init_func=init, frames=frames+1, interval=20)



    


















    
Out[8]:











    
    

    
    

    
    
    
    
    
    
    
    
    
  

     Once 
     Loop 
     Reflect 
  























In [ ]:



    






    











In [8]:



    


from sympy import symbols,Function,diff,exp



    











In [10]:



    


x,t,C =symbols("x t C")
f=Function("f")
c=5/sqrt(6)



    











In [12]:



    


N=1/(1+C*exp(-(x+c*t)/sqrt(6)))/(1+C*exp(-(x+c*t)/sqrt(6)))



    











In [14]:



    


leftT=diff(N,t);leftT



    


















    
Out[14]:








1.66666666666667*C*exp(-0.833333333333333*t - 0.408248290463863*x)/(C*exp(-0.833333333333333*t - 0.408248290463863*x) + 1)**3

















In [16]:



    


rightR=r*N*(1-N)+diff(N,x,2);rightR



    


















    
Out[16]:








C*(1.0*C*exp(-1.66666666666667*t - 0.816496580927726*x)/(C*exp(-0.833333333333333*t - 0.408248290463863*x) + 1) - 0.333333333333333*exp(-0.833333333333333*t - 0.408248290463863*x))/(C*exp(-0.833333333333333*t - 0.408248290463863*x) + 1)**3

















In [ ]:
$$u_t = d u_{xx} + B u (1 - u)$$ 




In [6]:



    


from scipy.integrate import simps
from scipy.sparse import spdiags, coo_matrix, csc_matrix, dia_matrix, dok_matrix, identity
from scipy.sparse.linalg import spsolve
from itertools import chain
import pylab 
import scipy



    











In [7]:



    


def M(x):
    """Returns coefficient matrix for diffusion equation (constant coefficients)"""
    ones = np.ones(len(x))
    M =  spdiags( [ones, -2*ones, ones], (-1,0,1), len(x), len(x) ).tocsc()
    J = len(x)
    M[J-1,J-2] = 2.0 # Neumann conditions
    return M
    
def set_DBCs(M):
    M = M.todok()
    M[0,0] = 1;     M[0,1] = 0;
    M[1,0] = 0;
    M[-2,-1] = 0
    M[-1,-1] = 1;   M[-1,-2] = 0
    return M.tocsc()



    











In [8]:



    


def bc_vector(x, d, h, left_value, right_value):
    bc = dok_matrix((len(x), 1)) # column vector
    bc[0,0]  = -left_value
    bc[1,0]  = d/h**2 * left_value
    bc[-2,0] = d/h**2 * right_value
    bc[-1,0] = -right_value
    return bc.tocsc()
 
def Imod(x):
    data = np.ones(len(x))
    data[0] = 0
    data[-1] = 0
    offset = 0
    shape = (len(x),len(x))
    Imod = dia_matrix((data, offset), shape=shape).todok()
    return Imod.tocsc()



    











In [9]:



    


def solution_fishers_eqn(x,t,d,B,n):
    """Analytical solution to Fishers equation with u(x0,t)=1
    and u_x(xJ,t)=0. The parameter n can be either 1 or 2."""
    if n == 1:
        U = 5 * np.sqrt(d * B / 6)
        return 1.0 / (1 + np.exp(np.sqrt(B/(6*d))*(x - U*t)) )**2
    elif n == 2:
        V = np.sqrt(d*B/2)
        return 1.0/ (1 + np.exp(np.sqrt(B/(2*d))*(x - V*t)) )
    else:
        raise ValueError("The parameter `n` must be an integer with the value 1 or 2.")



    











In [10]:



    


# Coefficients and simulation constants
x = np.linspace(-10, 10, 100)
x_range = np.max(x) - np.min(x)
h = x_range / len(x)
d = 0.1
B = 1.0
s = d / h**2
tau = 0.5
fsn = 1
 
# Initial conditions (needs to end up as a column vector)
u_init = np.matrix( solution_fishers_eqn(x,0.0,d,B,fsn) ).reshape((len(x),1))
u = u_init



    











In [11]:



    


# Matrices and vectors
MM = M(x)                # Coefficient matrix
Im = Imod(x)            # Modifed identity matrix
I = identity(len(x))    # Standard identity matrix
theta = 0.5 * np.ones(len(x))   # IMEX but with a fully
theta[0]  = 0.0                   # implicit at boundary
theta[-1] = 0.0                   # points.
theta_matrix = dia_matrix((theta, 0), shape=(len(x),len(x))).tocsc()    
theta_matrix_1m = dia_matrix((1-theta, 0), shape=(len(x),len(x))).tocsc()
 
# Boundary conditions
#M = set_DBCs(M)
left_value = 1.0
right_value = 0.0
bc = bc_vector(x, d, h, left_value, right_value)
print "left_value", left_value
print "right_value", right_value



    


















    






left_value 1.0
right_value 0.0

















In [12]:



    


# Cache
slices = [x]
 
# Time stepping loop
STEPS = 20
PLOTS = 5
for i in range(0, STEPS+1):
    
    u_array = np.asarray(u).flatten()
    # Reaction (nonlinear, calculated everytime)
    r = np.matrix(B*u_array**(fsn)*(1 - u_array)).reshape((len(x),1))
    
    # Export data
    
    if i % abs(STEPS/PLOTS) == 0:
        
        plot_num = abs(i/(STEPS/PLOTS))
        fraction_complete = plot_num / PLOTS
        print "fraction complete", fraction_complete
        
        print "I: %g" % (simps(u_array,dx=h), )
        pylab.figure(1)
        pylab.plot(x,u_array, "o", color=pylab.cm.Accent(fraction_complete))
        pylab.plot(x, solution_fishers_eqn(x, i*tau, d, B, fsn), "-", color=pylab.cm.Accent(fraction_complete))
        slices.append(u_array)
    
    # Prepare l.h.s. @ t_{n+1}
    A = (I - tau*s*theta_matrix*MM)
    
    # Prepare r.h.s. @ t_{n}
    b = csc_matrix( (I + tau*s*theta_matrix_1m*MM)*u + tau*r)
    
    # Set boundary conditions
    b[0,0] = left_value
    #b[-1,-1] = right_value
    
    # Solve linear 
    u = spsolve(A,b)                       # Returns an numpy array,
    u = np.matrix(u).reshape((len(x),1))   # need a column matrix



    


















    






fraction complete 0
I: 9.13315
fraction complete 0
I: 10.3442
fraction complete 0
I: 11.452
fraction complete 0
I: 12.5307
fraction complete 0
I: 13.6028
fraction complete 1
I: 14.6742










    
























In [ ]:

Gray-Scott Models:

\begin{eqnarray*} \frac{\partial U}{\partial t} & = & r_u\nabla^2 U -u v^2 +f(1-u) \cr \frac{\partial V}{\partial t} & = & r_v\nabla^2 V +u v^2 -(f+k)v \end{eqnarray*}

where $U, V$ is chemical spcies
Ref: 

1. Complex Patterns in a Simple System, John E. Pearson, Science 261, 5118, 189-192, 1993.
2.  Gray Scott Model of Reaction Diffusion, Abelson, Adams, Coore, Hanson, Nagpal, Sussman, http://http://groups.csail.mit.edu/mac/projects/amorphous/GrayScott/
3. Nicolas P. Rougier, http://www.loria.fr/~rougier/teaching/index.html




In [3]:



    


import numpy as np
from math import sin,cos,pi 
from scipy.integrate import ode,odeint
import matplotlib.pylab as plt
from matplotlib import animation
from JSAnimation import IPython_display
from ipywidgets import StaticInteract, RangeWidget,RadioWidget
from IPython.display import clear_output
import time,random
from scitools.all import movie
from IPython.html import widgets # Widget definitions
from IPython.display import display # Used to display widgets in the notebook
import os



    











In [5]:



    


directory = "imgs"
if not os.path.exists(directory):
    os.makedirs(directory)



    











In [6]:



    


n  = 200
#Du, Dv, F, k = 0.16, 0.08, 0.035, 0.065 # Bacteria 1
# Du, Dv, F, k = 0.14, 0.06, 0.035, 0.065 # Bacteria 2
# Du, Dv, F, k = 0.16, 0.08, 0.060, 0.062 # Coral
# Du, Dv, F, k = 0.19, 0.05, 0.060, 0.062 # Fingerprint
# Du, Dv, F, k = 0.10, 0.10, 0.018, 0.050 # Spirals
# Du, Dv, F, k = 0.12, 0.08, 0.020, 0.050 # Spirals Dense
# Du, Dv, F, k = 0.10, 0.16, 0.020, 0.050 # Spirals Fast
# Du, Dv, F, k = 0.16, 0.08, 0.020, 0.055 # Unstable
# Du, Dv, F, k = 0.16, 0.08, 0.050, 0.065 # Worms 1
# Du, Dv, F, k = 0.16, 0.08, 0.054, 0.063 # Worms 2
#Du, Dv, F, k = 0.16, 0.08, 0.035, 0.060 # Zebrafish
parameters=np.array([0.16, 0.08, 0.050, 0.065])



    











In [7]:



    


def GrayScott(parameters,title):
    
    Du=parameters[0];Dv=parameters[1];F=parameters[2];k=parameters[3];
    Z = np.zeros((n+2,n+2), [('U', np.double), ('V', np.double)])
    U,V = Z['U'], Z['V']
    u,v = U[1:-1,1:-1], V[1:-1,1:-1]

    r = 20
    u[...] = 1.0
    U[n/2-r:n/2+r,n/2-r:n/2+r] = 0.50
    V[n/2-r:n/2+r,n/2-r:n/2+r] = 0.25
    u += 0.05*np.random.random((n,n))
    v += 0.05*np.random.random((n,n))
    
    size = np.array(Z.shape)
    dpi = 72.0
    figsize= size[1]/float(dpi),size[0]/float(dpi)
    
    fig = plt.figure(figsize=figsize, dpi=dpi, facecolor="white")
    fig.suptitle(title,color='red',fontsize=16)
    fig.add_axes([0.0, 0.0, 1.0, 1.0], frameon=False)
    #im = plt.imshow(V, interpolation='bicubic', cmap=plt.cm.gray_r)
    im = plt.imshow(V, interpolation='bicubic', cmap=plt.cm.BuGn)
    plt.xticks([]), plt.yticks([])
    
    for i in xrange(25000):
        Lu = (U[0:-2,1:-1] +U[1:-1,0:-2] - 4*U[1:-1,1:-1] + U[1:-1,2:] +U[2:  ,1:-1] )
        Lv = (V[0:-2,1:-1] +V[1:-1,0:-2] - 4*V[1:-1,1:-1] + V[1:-1,2:] +V[2:  ,1:-1] )

        uvv = u*v*v
        u += (Du*Lu - uvv +  F *(1-u))
        v += (Dv*Lv + uvv - (F+k)*v)

        if i % 100 == 0:
           im.set_data(V)
           im.set_clim(vmin=V.min(), vmax=V.max())
           #plt.draw()
           # To make movie
           plt.savefig("imgs/tmp-%03d.png" % (i/100) ,dpi=dpi)
    movie('imgs/tmp*.png',encoder='ffmpeg', fps=20, vcodec='libvpx', quiet=True,
           output_file='GrayScott.webm')



    











In [8]:



    


GrayScott(parameters,title='')



    


















    








Found 250 files of the format imgs/tmp*.png.










    
























In [9]:



    


from IPython.core.display import HTML
video = open("GrayScott.webm", "rb").read()
video_encoded = video.encode("base64")
video_tag = '<video controls alt="test" src="data:video/webm;base64,{0}">'.format(video_encoded)
HTML(data=video_tag)



    


















    
Out[9]:






















In [46]:



    






    











In [10]:



    


mysecondwidget = widgets.RadioButtonsWidget(values=["Bacteria1","Bacteria2", "Coral","Fingerprint",  
                                            "Spirals", "Sprials Dense","Sprials Fast","Unstable",
                                            "Worms1", "Worms2","Zebrafish"],value="Worms1")
display(mysecondwidget)



    











In [17]:



    


if (mysecondwidget.value=='Bacteria1'):
     parameter= np.array([0.16, 0.08, 0.035, 0.065]);
elif (mysecondwidget.value=='Bacteria2'):
     parameter= np.array([0.14, 0.06, 0.035, 0.065]);        
elif (mysecondwidget.value=='Coral'):
     parameters= np.array([0.16, 0.08, 0.054, 0.063 ]);         
elif (mysecondwidget.value=='Fingerprint'):
     parameters= np.array([0.19, 0.05, 0.060, 0.062  ]); 
elif (mysecondwidget.value=='Sprials'):
     parameters= np.array([0.10, 0.10, 0.018, 0.050  ]);
elif (mysecondwidget.value=='Sprials Dense'):
     parameters= np.array([0.12, 0.08, 0.020, 0.050 ]); 
elif (mysecondwidget.value=='Sprials Fast'):
     parameters= np.array([0.10, 0.16, 0.020, 0.050 ]);  
elif (mysecondwidget.value=='Unstable'):
     parameters= np.array([0.16, 0.08, 0.020, 0.055  ]); 
elif (mysecondwidget.value=='Worms1'):
     parameter= np.array([0.16, 0.08, 0.050, 0.065]);
elif (mysecondwidget.value=='Worms2'):
     parameters= np.array([0.16, 0.08, 0.054, 0.063 ]);        
elif (mysecondwidget.value=='Zebrafish'):
     parameters= np.array([0.16, 0.08, 0.035, 0.060 ]);



    











In [18]:



    


GrayScott(parameters,mysecondwidget.value)



    


















    








Found 250 files of the format imgs/tmp*.png.










    
























In [19]:



    


from IPython.core.display import HTML
video = open("GrayScott.webm", "rb").read()
video_encoded = video.encode("base64")
video_tag = '<video controls alt="test" src="data:video/webm;base64,{0}">'.format(video_encoded)
HTML(data=video_tag)



    


















    
Out[19]:






















In [121]:



    






    











In [178]:



    


fwidget = widgets.FloatSliderWidget(description='f value:', min='0.01', max='0.05', step='0.001', value='0.024')
kwidget = widgets.FloatSliderWidget(description='k value:', min='0.04', max='0.07', step='0.005', value='0.055')



    











In [179]:



    


display(fwidget)
display(kwidget)



    











In [180]:



    


parameters=np.array([fwidget.value,kwidget.value])



    











In [181]:



    


parameters



    


















    
Out[181]:








array([ 0.035,  0.065])

















In [187]:



    


def GrayScott2(parameters):
    
    Du=0.01;Dv=0.005;F=parameters[0];k=parameters[1];
    Z = np.zeros((n+2,n+2), [('U', np.double), ('V', np.double)])
    U,V = Z['U'], Z['V']
    u,v = U[1:-1,1:-1], V[1:-1,1:-1]

    r = 4
    u[...] = 1.0
    U[n/2-r:n/2+r,n/2-r:n/2+r] = 0.50
    V[n/2-r:n/2+r,n/2-r:n/2+r] = 0.25
    #u += 0.05*np.random.random((n,n))
    #v += 0.05*np.random.random((n,n))
    
    size = np.array(Z.shape)
    dpi = 72.0
    figsize= size[1]/float(dpi),size[0]/float(dpi)
    
    fig = plt.figure(figsize=figsize, dpi=dpi, facecolor="white")
    fig.add_axes([0.0, 0.0, 1.0, 1.0], frameon=False)
    #im = plt.imshow(V, interpolation='bicubic', cmap=plt.cm.gray_r)
    im = plt.imshow(V, interpolation='bicubic', cmap=plt.cm.BuGn)
    plt.xticks([]), plt.yticks([])
    
    for i in xrange(25000):
        Lu = (U[0:-2,1:-1] +U[1:-1,0:-2] - 4*U[1:-1,1:-1] + U[1:-1,2:] +U[2:  ,1:-1] )
        Lv = (V[0:-2,1:-1] +V[1:-1,0:-2] - 4*V[1:-1,1:-1] + V[1:-1,2:] +V[2:  ,1:-1] )

        uvv = u*v*v
        u += (Du*Lu - uvv +  F *(1-u))
        v += (Dv*Lv + uvv - (F+k)*v)

        if i % 100 == 0:
           im.set_data(V)
           im.set_clim(vmin=V.min(), vmax=V.max())
           #plt.draw()
           # To make movie
           plt.savefig("imgs/tmp-%03d.png" % (i/100) ,dpi=dpi)
    movie('imgs/tmp*.png',encoder='ffmpeg', fps=20, vcodec='libvpx', quiet=True,
           output_file='GrayScott.webm')



    











In [188]:



    


GrayScott2(parameters)



    


















    








Found 250 files of the format imgs/tmp*.png.










    
























In [184]:



    


from IPython.core.display import HTML
video = open("GrayScott.webm", "rb").read()
video_encoded = video.encode("base64")
video_tag = '<video controls alt="test" src="data:video/webm;base64,{0}">'.format(video_encoded)
HTML(data=video_tag)



    


















    
Out[184]:






















In [65]:



    


# Convert to HTML format
# Jake VanderPlas
# Director of Research – Physical Sciences
# eScience Institute, University of Washington

from IPython.nbformat import current as nbformat
from IPython.nbconvert.exporters import HTMLExporter

input_file = 'diffusion.ipynb'
output_file = 'diffusion.html'

with open(input_file) as f:
    notebook_content = f.read()

exporter = HTMLExporter(template_file='full')
nb_json = nbformat.reads_json(notebook_content)
(body, resources) = exporter.from_notebook_node(nb_json)

with open(output_file, 'w') as f:
    f.write(body)



    


















    






/Users/cch/anaconda/lib/python2.7/site-packages/IPython/nbconvert/filters/markdown.py:78: UserWarning: Node.js 0.9.12 or later wasn't found.
Nbconvert will try to use Pandoc instead.
  "Nbconvert will try to use Pandoc instead.")

















In [ ]:



    


%%bash

# upload Python source
cd ~
# extract source
cd src
sage --python setup.py install --user

Content source: cchuang2009/DyLabIPy

Similar notebooks:

notebook.community | gallery | about