3D (Polar/Cylindrical Coordinate) Animation of 2D Diffusion Equation using Python, Scipy, and Matplotlib

 Yup, that same code but in polar coordinate.

 I use nabla operator for cylindrical coordinate but ditch the z component.

 So, what's the z-axis for? It's represent the u value, in this case, temperature, as function of r and phi (I know I should use rho, but, ...)

import scipy as sp

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import matplotlib.pyplot as plt

import mpl_toolkits.mplot3d.axes3d as p3
import matplotlib.animation as animation

#dr  = .1
#dp  = .1
#nr      = int(1/dr)
#np      = int(2*sp.pi/dp)


nr  = 10
np  = 10

r       = sp.linspace(0.,1.,nr)
p       = sp.linspace(0.,2*sp.pi,np)

dr  = r[1]-r[0]
dp  = p[1]-p[0]

a   = .5

tmax    = 100
t       = 0.


dr2     = dr**2
dp2     = dp**2

dt      = dr2 * dp2 / (2 * a * (dr2 + dp2) )
dt      /=10.
print 'dr = ',dr 
print 'dp = ',dp 
print 'dt = ',dt



ut      = sp.zeros([nr,np])
u0      = sp.zeros([nr,np])

ur      = sp.zeros([nr,np])
ur2     = sp.zeros([nr,np])


#initial

for i in range(nr):
    for j in range(np):
        if ((i>4)&(i<6)):
            u0[i,j] = 1.
#print u0

def hitung_ut(ut,u0):
    for i in sp.arange (len(r)):
        if r[i]!= 0.:
            ur[i,:]     = u0[i,:]/r[i]
            ur2[i,:]     = u0[i,:]/(r[i]**2)
    ut[1:-1, 1:-1]  = u0[1:-1, 1:-1] + a*dt*(
            (ur[1:-1, 1:-1] - ur[:-2, 1:-1])/dr+
            (u0[2:, 1:-1] - 2*u0[1:-1, 1:-1] + u0[:-2,1:-1])/dr2+
            (ur2[1:-1, 2:] - 2*ur2[1:-1, 1:-1] + ur2[1:-1, :-2])/dp2)
    #calculate the edge
    ut[1:-1, 0]  = u0[1:-1, 0] + a*dt*(
            (ur[1:-1, 0] - ur[:-2, 0])/dr+
            (u0[2:, 0] - 2*u0[1:-1, 0] + u0[:-2, 0])/dr2+
            (ur2[1:-1, 1] - 2*ur2[1:-1, 0] + ur2[1:-1, np-1])/dp2)
    ut[1:-1, np-1]  = u0[1:-1, np-1] + a*dt*(
            (ur[1:-1, np-1] - ur[:-2, np-1])/dr+
            (u0[2:, np-1] - 2*u0[1:-1, np-1] + u0[:-2,np-1])/dr2+
            (ur2[1:-1, 0] - 2*ur2[1:-1, np-1] + ur2[1:-1, np-2])/dp2)


#hitung_ut(ut,u0)
#print ut

def data_gen(framenumber, Z ,surf):
    global ut
    global u0
    global t
    hitung_ut(ut,u0)
    u0[:] = ut[:]
    Z = u0
    t += 1
    print t
    
    ax.clear()
    plotset()
    surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,
                       linewidth=0, antialiased=False, alpha=0.7)
    return surf,


fig = plt.figure()
#ax = fig.gca(projection='3d')
ax = fig.add_subplot(111, projection='3d')

P,R = sp.meshgrid(p,r)

X,Y = R*sp.cos(P),R*sp.sin(P) 

Z = u0
print len(R), len(P)



def plotset():
    ax.set_xlim3d(-1., 1.)
    ax.set_ylim3d(-1., 1.)
    ax.set_zlim3d(-1.,1.)
    ax.set_autoscalez_on(False)
    ax.zaxis.set_major_locator(LinearLocator(10))
    ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
    cset = ax.contour(X, Y, Z, zdir='x', offset=-1. , cmap=cm.coolwarm)
    cset = ax.contour(X, Y, Z, zdir='y', offset=1. , cmap=cm.coolwarm)
    cset = ax.contour(X, Y, Z, zdir='z', offset=-1., cmap=cm.coolwarm)

plotset()
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,
                       linewidth=0, antialiased=False, alpha=0.7)

fig.colorbar(surf, shrink=0.5, aspect=5)

ani = animation.FuncAnimation(fig, data_gen, fargs=(Z, surf),frames=4096, interval=4, blit=False)
#ani.save('2dDiffusionfRadialf1024b512.mp4', bitrate=1024)

plt.show()    

.

100x100 size

The Wrong Code Will often Provide Beautiful Result, :)

 It means to compute 2d diffusion equation just like previous post in polar/cylindrical coordinate, and all went to wrong direction, :)

 Still trying to understand matplotlib mplot3d behavior

import scipy as sp

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import matplotlib.pyplot as plt

import mpl_toolkits.mplot3d.axes3d as p3
import matplotlib.animation as animation

#dr  = .1
#dp  = .1
#nr      = int(1/dr)
#np      = int(2*sp.pi/dp)

nr  = 10
np  = 10

dr  = 1./nr
dp  = 2*sp.pi/np

a   = .5

tmax    = 100
t       = 0.


dr2     = dr**2
dp2     = dp**2

dt      = dr2 * dp2 / (2 * a * (dr2 + dp2) )
dt      /=10.
print dt



ut      = sp.zeros([nr,np])
u0      = sp.zeros([nr,np])

ur      = sp.zeros([nr,np])
ur2     = sp.zeros([nr,np])

r       = sp.arange(0.,1.,dr)
p       = sp.arange(0.,2*sp.pi,dp)

#initial

for i in range(nr):
    for j in range(np):
        if ( (i>(2*nr/5.)) & (i<(3.*nr/3.)) ):
            u0[i,j] = 1.
#print u0

def hitung_ut(ut,u0):
    for i in sp.arange (len(r)):
        if r[i]!= 0.:
            ur[i,:]     = u0[i,:]/r[i]
            ur2[i,:]     = u0[i,:]/(r[i]**2)
    ut[1:-1, 1:-1]  = u0[1:-1, 1:-1] + a*dt*(
            (ur[1:-1, 1:-1] - ur[:-2, 1:-1])/dr+
            (u0[2:, 1:-1] - 2*u0[1:-1, 1:-1] + u0[:-2,1:-1])/dr2+
            (ur2[1:-1, 2:] - 2*ur2[1:-1, 1:-1] + ur2[1:-1, :-2])/dp2)


#hitung_ut(ut,u0)
#print ut


def data_gen(framenumber, Z ,surf):
    global ut
    global u0
    hitung_ut(ut,u0)
    u0[:] = ut[:]
    Z = u0
    
    ax.clear()
    plotset()
    surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,
                       linewidth=0, antialiased=False, alpha=0.7)
    return surf,


fig = plt.figure()
#ax = fig.gca(projection='3d')
ax = fig.add_subplot(111, projection='3d')

R = sp.arange(0,1,dr)
P = sp.arange(0,2*sp.pi,dp)
R,P = sp.meshgrid(R,P)

X,Y = R*sp.cos(P),R*sp.sin(P) 

Z = u0
print len(R), len(P)



def plotset():
    ax.set_xlim3d(-1., 1.)
    ax.set_ylim3d(-1., 1.)
    ax.set_zlim3d(-1.,1.)
    ax.set_autoscalez_on(False)
    ax.zaxis.set_major_locator(LinearLocator(10))
    ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
    cset = ax.contour(X, Y, Z, zdir='x', offset=0. , cmap=cm.coolwarm)
    cset = ax.contour(X, Y, Z, zdir='y', offset=1. , cmap=cm.coolwarm)
    cset = ax.contour(X, Y, Z, zdir='z', offset=-1., cmap=cm.coolwarm)

plotset()
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,
                       linewidth=0, antialiased=False, alpha=0.7)

fig.colorbar(surf, shrink=0.5, aspect=5)

ani = animation.FuncAnimation(fig, data_gen, fargs=(Z, surf),frames=500, interval=30, blit=False)
#ani.save('2dDiffusionf500b512.mp4', bitrate=512)

plt.show()    

.

Miss Lay.

 Baca komen seperti  "Permainan gitar Eross pecah banget...." di Youtube, hati langsung panas.

 Berani-beraninya menghina sang maestro, sana, ke comberan sana, tempat asalmu, ...

 Eh, tunggu...

 Itu bukan hinaan, itu pujian...

 Pecah? Pujian?

 Yeah, beberapa dari kita sering menggunakan kata tak pada tempatnya, :)

 Penggunaan kata 'pecah' sebagai pujian datang dari Maia Estianty saat menjadi juri The Remix.

 Tentu saja penggunaan kata 'pecah' di dunia DJ atau remix sangat cocok; merujuk pada sunyi yang lambat laun menjadi semakin ramai seiring dengan penumpukan nada-nada sampling yang makin banyak, tempo yang makin cepat dan diakhiri dengan sebuah hentakan..., pecah..., 

 OK, itu pujian, bagaimana kalo kata 'pecah' digunakan di permainan gitar. Yeah, kami menganggap kata pecah sebagai kata berkonotasi buruk; merujuk pada suara gitar yang seharusnya clean namun karena setting yang buruk jadi terdengar 'brebet' atau 'pecah'.

 Bukannya efek gitar Telecaster-nya Eross mempunya karakter pecah?  Ehm, EDrive-nya Eross, seperti namanya, 'drive', hanya bertugas sebagai booster saat permainan lead, bahkan kalo kita lihat di penampilan live, Eross jarang sekali menggunakan stompbox, dia memainkan knob volume untuk efek crunch atau crisp. Jika  ingin suara clean, dia mengecilkan volume gitarnya.

 Suara Telecaster Eross 'pecah'? Well, lebih baik pakai kata crunchy, crispy, ..., lebih diterima oleh banyak gitaris, :)

3D Animation of 2D Diffusion Equation using Python, Scipy, and Matplotlib

 I wrote the code on OS X El Capitan, use a small mesh-grid.  Basically it's same code like the previous post.

 I use surface plot mode for the graphic output and animate it.

 Because my Macbook Air is suffered from running laborious code, I save the animation on my Linux environment, 1024 bitrate, 1000 frames.

story

import scipy as sp
import time
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import matplotlib.pyplot as plt

import mpl_toolkits.mplot3d.axes3d as p3
import matplotlib.animation as animation



dx=0.01        
dy=0.01     
a=0.5          
timesteps=500  
t=0.

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


dx2=dx**2 
dy2=dy**2 

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

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

for i in range(nx):
 for j in range(ny):
  if ( ( (i*dx-0.5)**2+(j*dy-0.5)**2 <= 0.1)
   & ((i*dx-0.5)**2+(j*dy-0.5)**2>=.05) ):
    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 )
        

def data_gen(framenumber, Z ,surf):
    global u
    global ui
    evolve_ts(u,ui)
    ui[:] = u[:]
    Z = ui
    
    ax.clear()
    plotset()
    surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,
                       linewidth=0, antialiased=False, alpha=0.7)
    return surf,


fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

X = sp.arange(0,1,dx)
Y = sp.arange(0,1,dy)
X,Y= sp.meshgrid(X,Y)

Z = ui 

def plotset():
    ax.set_xlim3d(0., 1.)
    ax.set_ylim3d(0., 1.)
    ax.set_zlim3d(-1.,1.)
    ax.set_autoscalez_on(False)
    ax.zaxis.set_major_locator(LinearLocator(10))
    ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
    cset = ax.contour(X, Y, Z, zdir='x', offset=0. , cmap=cm.coolwarm)
    cset = ax.contour(X, Y, Z, zdir='y', offset=1. , cmap=cm.coolwarm)
    cset = ax.contour(X, Y, Z, zdir='z', offset=-1., cmap=cm.coolwarm)

plotset()
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,
                       linewidth=0, antialiased=False, alpha=0.7)


fig.colorbar(surf, shrink=0.5, aspect=5)

ani = animation.FuncAnimation(fig, data_gen, fargs=(Z, surf),frames=1000, interval=30, blit=False)
ani.save("2dDiffusion.mp4", bitrate=1024)

#plt.show()    

.

2D Diffusion Equation using Python, Scipy, and VPython

 I got it from here, but modify it here and there.

 I also add animation using vpython but can't find 3d or surface version, so I planned to go to matplotlib surface plot route, :)

 (update: here it is, :) )

#!/usr/bin/env python
"""
A program which uses an explicit finite difference
scheme to solve the diffusion equation with fixed
boundary values and a given initial value for the
density.

Two steps of the solution are stored: the current
solution, u, and the previous step, ui. At each time-
step, u is calculated from ui. u is moved to ui at the
end of each time-step to move forward in time.

http://www.timteatro.net/2010/10/29/performance-python-solving-the-2d-diffusion-equation-with-numpy/

he uses matplotlib

I use visual python

"""
import scipy as sp
import time
from visual import *
from visual.graph import *


graph1 = gdisplay(x=0, y=0, width=600, height=400, 
          title='x vs. T', xtitle='x', ytitle='T', 
          foreground=color.black, background=color.white)


# Declare some variables:

dx=0.01        # Interval size in x-direction.
dy=0.01        # Interval size in y-direction.
a=0.5          # Diffusion constant.
timesteps=500  # Number of time-steps to evolve system.
t=0.

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

dx2=dx**2 # To save CPU cycles, we'll compute Delta x^2
dy2=dy**2 # and Delta y^2 only once and store them.

# For stability, this is the largest interval possible
# for the size of the time-step:
dt = dx2*dy2/( 2*a*(dx2+dy2) )

# Start u and ui off as zero matrices:
ui = sp.zeros([nx,ny])
u = sp.zeros([nx,ny])

# Now, set the initial conditions (ui).
for i in range(nx):
 for j in range(ny):
  if ( ( (i*dx-0.5)**2+(j*dy-0.5)**2 <= 0.1)
   & ((i*dx-0.5)**2+(j*dy-0.5)**2>=.05) ):
    ui[i,j] = 1
'''
def evolve_ts(u, ui):
 global nx, ny
 """
 This function uses two plain Python loops to
 evaluate the derivatives in the Laplacian, and
 calculates u[i,j] based on ui[i,j].
 """
 for i in range(1,nx-1):
  for j in range(1,ny-1):
   uxx = ( ui[i+1,j] - 2*ui[i,j] + ui[i-1, j] )/dx2
   uyy = ( ui[i,j+1] - 2*ui[i,j] + ui[i, j-1] )/dy2
   u[i,j] = ui[i,j]+dt*a*(uxx+uyy)
'''
def evolve_ts(u, ui):
 """
 This function uses a numpy expression to
 evaluate the derivatives in the Laplacian, and
 calculates u[i,j] based on ui[i,j].
 """
 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 )
        
# Now, start the time evolution calculation...
#tstart = time.time()

f1 = gcurve(color=color.blue)

while True:
    rate(60)
    #for m in range(1, timesteps+1):
    if t<timesteps:
        t+=dt
 evolve_ts(u, ui)
        ui[:] = u[:] # I add this line to update ui value (not present in original code)
 #print "Computing u for m =", m
    f1.gcurve.pos   =   []
    for i in arange(nx):
        f1.plot(pos=(i,u[nx/2,i]))
    
#tfinish = time.time()

#print "Done."
#print "Total time: ", tfinish-tstart, "s"
#print "Average time per time-step using numpy: ", ( tfinish - tstart )/timesteps, "s."

.

Numpy Slice Expression

 Suppossed we have two array a and b

 If we want to set b as finite difference result of a, we may tempted to do this

for i in range (9):
 b[i] = a[i+1]-a[i]

There's another (faster) way. The performance's close to the pure C, :)

b[:-1] = a[1:]-a[:-1]

What's that?

Numpy has slice form for array. If we have an array with length 10, the a[:] refers to all value in a.

a[1:] refers to a[1] to a[9] (without a[0])
a[3:] refers to a[3] to a[9]
a[:-1] refers to a[0] to a[8]
a[:-3] refers to a[0] to a[6]
a[1:-1] refers to a[1] to a[8]
...
and so on

Here's my tinkering with slice expression

>>> from numpy import *
>>> a = zeros(10)
>>> b = zeros(10)
>>> a[5]=1.
>>> a
array([ 0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.])
>>> b
array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.])
>>> a[6]=2.
>>> a
array([ 0.,  0.,  0.,  0.,  0.,  1.,  2.,  0.,  0.,  0.])
>>> b[:-1]=a[:-1]-a[1:]
>>> b
array([ 0.,  0.,  0.,  0., -1., -1.,  2.,  0.,  0.,  0.])
>>> b[:-1]=a[:-1]+a[1:]
>>> b
array([ 0.,  0.,  0.,  0.,  1.,  3.,  2.,  0.,  0.,  0.])
>>> 

I like Python, :)

Car Free Day

 Sepertinya ada (banyak) yang mengartikan sebagai kendaraan bebas berjalan di manapun, bahkan di tempat yang saat normal gak bisa dimasuki, :(

 (hampir ketabrak sepeda motor yang dikendarai mahasiswa berboncengan  di jalan antara gedung MIPA dan Fisika)

 #EdisiError