Aravir

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.

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() .

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.

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.

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.