Python画图

Python画图

全栈程序员-用户IM • 2022年7月5日下午11:36 • 未分类

大家好，又见面了，我是你们的朋友全栈君。

帮别人实现rubber friction的运算公式。公式不是我专业的内容，完全不懂，只是照着公式用python计算并画图。做出的图，与论文中的图进行对比，完全一致。

之所以用python，是因为matlab实在太大了。。。

python科学计算需要用到numpy和scipy，画图需要用到matplotlib。这三个模块的安装，windows下很容易，mac下稍微麻烦一些。

matplotlib需要注意的是，因为公式涉及的三个图，分别是对数-对数，对数-常数，常数-对数坐标，所以不能仅仅使用plot。

步骤一：

公式： Python画图

代码：

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)

图：

Python画图

步骤二：

公式： Python画图（注：公式有些错误，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)