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 |
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()
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()
Add a comment