更新空いてしまいました。
Parabolic Bowl
先週末から平面二次元の計算を始めました。 手始めにparabolic bowlの図化をpythonで書いてみました。 気が向いたら回してみてください。DBに置いておきます。 さすがに遅すぎるので実計算には使えそうにないです。動画っぽくしたいなら、ImageMagickでgifにするほうが良いかもしれないです。
ソースコードはこんな感じです。ほぼグラフのところで、解析解の計算は数行です。やっぱりI/Oは手間ですね。
# -*- coding: UTF-8 -*- # use python3.4 ## @package # # # ###%matplotlib inline from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt from matplotlib import cm import numpy as np class MyGragh: def __init__(self,X,Y,zb,wL,idX,idY): self.X,self.Y,self.idX,self.idY = X,Y,idX,idY #make gragh self.fig = plt.figure(figsize=(12, 10)) #make 3Dgragh self.ax = self.fig.add_subplot(2, 1, 1, projection='3d') self.ax.set_zlim3d(0, 14) self.ax.view_init(elev=33, azim=38) Xm, Ym = np.meshgrid(self.X, self.Y) h = list(map(list,zip(*zb))) self.ax.plot_surface(Xm, Ym, h, color='k', zorder=1, rstride=1, cstride=1,linewidth=0.1, antialiased=False,alpha=0.3) h = list(map(list,zip(*wL))) self.surf = self.ax.plot_surface(Xm, Ym, h, cmap=cm.jet, zorder=2, rstride=1, cstride=1,linewidth=0.1, antialiased=False, alpha=1) #make 2Dgragh self.ax2 = self.fig.add_subplot(2, 2, 3) plt.title('X = '+str(self.X[self.idX])) h = zb[self.idX,:] plt.plot(self.Y,h, color='k') h = wL[self.idX,:] self.lineX = self.ax2.plot(self.Y,h, color='b') self.ax3 = self.fig.add_subplot(2, 2, 4) plt.title('Y = '+str(self.Y[self.idY])) h = zb[:,self.idY] plt.plot(self.X,h, color='k') h = wL[:,self.idX] self.lineY = self.ax3.plot(self.X,h, color='b') def plotWL(self,wL,t,u,v): #make 3Dgragh self.surf.remove() Xm, Ym = np.meshgrid(self.X, self.Y) h = list(map(list,zip(*wL))) self.surf = self.ax.plot_surface(Xm, Ym, h, cmap=cm.jet, zorder=2, rstride=1, cstride=1,linewidth=0.1, antialiased=False, alpha=1) #make 2Dgragh for c in self.lineX : c.remove() h = wL[self.idX,] self.lineX = self.ax2.plot(self.Y,h, color='b') for c in self.lineY : c.remove() h = wL[:,self.idY] self.lineY = self.ax3.plot(self.X,h, color='b') #make title self.fig.suptitle('t='+str(t)+'[sec]'+' U='+str(round(u,3))+'[m/s]'+' V='+str(round(v,3))+'[m/s]') plt.pause(0.001) def __enter__(self): return self def __exit__(self, exc_type, exc_value, traceback): plt.savefig('output.png') return False # main nx, ny = 21, 21 X = np.arange(-10, 11, 1) Y = np.arange(-10, 11, 1) a = 7 h0 = 4 zb = np.array([[h0/a**2*(X[i]**2+Y[j]**2) for j in range(ny)] for i in range(nx)]) eta = a/3 tPriod=60 omega = 2*np.pi/tPriod #wL = lambda t: np.array([[ eta*h0/a**2*(2*X[i]*np.cos(t*omega)+2*Y[j]*np.sin(t*omega)-eta)+h0 for j in range(ny)] for i in range(nx)]) wL = lambda t, x, y: eta*h0/a**2*(2*x*np.cos(t*omega)+2*y*np.sin(t*omega)-eta)+h0 wLr = lambda t: np.array([[ wL(t,X[i],Y[j]) if wL(t,X[i],Y[j]) - zb[i,j] > 0 else zb[i,j] for j in range(ny)] for i in range(nx)]) u = lambda t: -eta*omega*np.sin(t*omega) v = lambda t: -eta*omega*np.cos(t*omega) idX=10 idY=10 tmax = 600 with MyGragh(X,Y,zb,wLr(0),idX,idY) as g: for t in range(1,tmax): if t%10==0 : g.plotWL(wLr(t),t,u(t),v(t))
アジョイント法
とりあえず、矩形断面ではできたっぽくて、テンションが上がってます。パッと見せれる絵がないので勉強会のときにでも話します。
