Python画图
帮别人实现rubber friction的运算公式。公式不是我专业的内容，完全不懂，只是照着公式用python计算并画图。做出的图，与论文中的图进行对比，完全一致。
之所以用python，是因为matlab实在太大了。。。
python科学计算需要用到numpy和scipy，画图需要用到matplotlib。这三个模块的安装，windows下很容易，mac下稍微麻烦一些。
matplotlib需要注意的是，因为公式涉及的三个图，分别是对数-对数，对数-常数，常数-对数坐标，所以不能仅仅使用plot。
步骤一：
公式：
代码：
View Code
1 def EOmega(E1,omega,tau,a1): 2 a=complex(1,-omega*tau) 3 b=complex(a1+1,-omega*tau) 4 return E1*a/b 5 6 def DrawEOmega(E1,tau,a1,xmin,xmax): 7 re=[] 8 im=[] 9 10 omega = np.logspace(xmin, xmax, 1000) 11 for o in omega: 12 re.append(EOmega(E1,o,tau,a1).real) 13 im.append(abs(EOmega(E1,o,tau,a1).imag)) 14 15 plt.figure() 16 plt.loglog(omega,re,label="$Re(E)$",color="red",linewidth=2) 17 plt.loglog(omega,im,"b--",label="$Im(E)$") 18 plt.ylim(10**0,10**10) 19 plt.xlabel("frequency omega[rad/s]") 20 plt.ylabel("E-modulus[Pa]") 21 plt.title("Complex E-modulus for rheological model") 22 plt.legend(loc="center right") 23 plt.show() 24 25 26 E1=10.0**9 27 a1=1000.0 28 tau=1.0/10**3 29 xmin=-2 30 xmax=8 31 DrawEOmega(E1,tau,a1,xmin,xmax)
图：
步骤二：
公式：
（注：公式有些错误，G(q)的第二个积分没有Im）
代码：
View Code
1 def EOmegaExpand(E1,tau,a1,q0,zeta,v,phi): 2 omega=q0*zeta*v*np.cos(phi) 3 return EOmega(E1,omega,tau,a1) 4 5 def GQPartTwo(phi,E1,tau,a1,q0,zeta,v,gamma,sigma0): 6 a=EOmegaExpand(E1,tau,a1,q0,zeta,v,phi) 7 b=(1.0-gamma**2)*sigma0 8 c=a/b 9 re=c.real 10 im=c.imag 11 d=re**2+im**2 12 return d 13 14 def GQPartOne(zeta,E1,tau,a1,q0,v,gamma,sigma0,H): 15 a=quad(GQPartTwo, 0, 2*np.pi,args=(E1,tau,a1,q0,zeta,v,gamma,sigma0))[0] 16 b=(zeta**(1.0-2*H))*a 17 return b 18 19 def GQ(E1,tau,a1,q0,v,sigma0,gamma,H,h0,q): 20 a=quad(GQPartOne,1,q/q0,args=(E1,tau,a1,q0,v,gamma,sigma0,H))[0] 21 b=((q0*h0)**2)*H 22 c=a*b/(16*np.pi) 23 return c 24 25 def PQ(E1,tau,a1,q0,sigma0,gamma,H,h0,v,zeta): 26 q=q0*zeta 27 a=GQ(E1,tau,a1,q0,v,sigma0,gamma,H,h0,q) 28 b=1/(2*np.sqrt(a)) 29 c=erf(b) 30 return c 31 32 def DrawPQ(E1,tau,a1,q0,sigma0,gamma,H,h0,xmin,xmax): 33 p1=[] 34 p2=[] 35 p3=[] 36 p4=[] 37 p5=[] 38 39 v1=0.00005 40 v2=0.013 41 v3=0.316 42 v4=7.85 43 v5=187 44 45 l = np.linspace(xmin, xmax, 50) 46 for o in l: 47 p1.append(PQ(E1,tau,a1,q0,sigma0,gamma,H,h0,v1,o)) 48 p2.append(PQ(E1,tau,a1,q0,sigma0,gamma,H,h0,v2,o)) 49 p3.append(PQ(E1,tau,a1,q0,sigma0,gamma,H,h0,v3,o)) 50 p4.append(PQ(E1,tau,a1,q0,sigma0,gamma,H,h0,v4,o)) 51 p5.append(PQ(E1,tau,a1,q0,sigma0,gamma,H,h0,v5,o)) 52 53 plt.figure() 54 plt.semilogy(l,p1,"b--",label="$v=0.00005m/s$",linewidth=5) 55 plt.semilogy(l,p2,label="$v=0.013m/s$",linewidth=5) 56 plt.semilogy(l,p3,"r--",label="$v=0.316m/s$",linewidth=1) 57 plt.semilogy(l,p4,"b--",label="$v=7.85m/s$",linewidth=1) 58 plt.semilogy(l,p5,label="$v=187m/s$",linewidth=1) 59 plt.xlabel("zeta[-]") 60 plt.ylabel("P[-]") 61 plt.title("Apparent contact area") 62 plt.legend(loc='center left', bbox_to_anchor=(1, 0.5)) 63 plt.show() 64 65 E1=10.0**9 66 a1=1000.0 67 tau=1.0/10**3 68 q0=2000.0 69 sigma0=0.2*(10**6) 70 gamma=0.5 71 H=0.85 72 h0=5.0/10**4 73 xmin=1 74 xmax=100 75 DrawPQ(E1,tau,a1,q0,sigma0,gamma,H,h0,xmin,xmax)
图：
步骤三：
公式：
这一步的代码没有什么难度，只是注意一下坐标的问题即可。
目前还有一点问题，即第二张图中的图例不知道为什么总是无法显示完全，需要再查一下官方文档。
