Python画图[通俗易懂]

Python画图[通俗易懂]帮别人实现rubberfriction的运算公式。公式不是我专业的内容,完全不懂,只是照着公式用python计算并画图。做出的图,与论文中的图进行对比,完全一致。之所以用python,是因为mat

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

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

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

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

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

步骤一:

公式:Python画图[通俗易懂]

代码:

Python画图[通俗易懂]
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)

代码:

Python画图[通俗易懂]
Python画图[通俗易懂]
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)

图:

Python画图[通俗易懂]

步骤三:

公式:Python画图[通俗易懂]

这一步的代码没有什么难度,只是注意一下坐标的问题即可。

 

目前还有一点问题,即第二张图中的图例不知道为什么总是无法显示完全,需要再查一下官方文档。

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。

发布者:全栈程序员-用户IM,转载请注明出处:https://javaforall.cn/155983.html原文链接:https://javaforall.cn

【正版授权,激活自己账号】: Jetbrains全家桶Ide使用,1年售后保障,每天仅需1毛

【官方授权 正版激活】: 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...

(0)
blank

相关推荐

  • urlopen error errno111(英雄联盟报错error)

    原因获取地址信息失败,通常是由于自动选择的DNS服务器不行解决方法更换DNS服务器控制面板(win+R+control)->网络和Internet2.打开网络共享中心3.打开当前连接的网络,以wifi为例4.打开属性5.打开Internet协议版本46.选择使用特定DNS服务器地址,并设置为114.114.114.1148.设置完成后,错误解决…

  • 大数据平台建设路径

    大数据平台建设路径数字化转型这几年成为了集团公司或大型企业寻求业务突破的重要手段,而且各行各业都不乏成功案例,例如金融行业的建设银行、招商银行、平安保险等,还有一些生动的例子,例如百丽鞋业。而国外的经典案例就是亚马逊和特斯拉。集团客户动静都很大,首先要做咨询规划,必须跟上投入和管理配套,然后可以开始平台建设,最终开始应用,赋能各业务线和产品线。从技术的角度,一个完整的大数据平台通常会三期建设。先看下图:蓝色部分是第一期内容,绿色部分是第二期内容,橙色部分是第三期内容。第一步,先建设大数据基础平台。首要满足的是数据汇聚

  • 根据sessionid获取session对象_sessionattributes注解

    根据sessionid获取session对象_sessionattributes注解session.setAttribute(“sessionName”,Object);用来设置session值的,sessionName是名称,object是你要保存的对象。session.getAttribute(“sessionName”);用来得到对应名称的session值,即得到object对象,注意需要进行类型转换!

  • 安全帽识别的原理和系统应用

    安全帽识别的原理和系统应用安全帽识别的原理是用AI技术对工作现场的视频进行实时分析,如果发现工作人员未按要求佩戴安全帽或违规吸烟,系统会自动发出警报,在提醒监理人员的同时,系统会自动保存时间、地点及相应的照片,作为处罚的依据,AI大潮之下,传统的建筑工地也有许多方面得到了提升,除了人脸识别之外,智慧工地最关心的就是安全问题,如何将AI技术应用于安全方面呢?鹰眸安全帽识别系统就是这些特殊区域的守护者。可以说,安全帽智能识别是…

  • Fiddler HPPTS 配置教程「建议收藏」

    Fiddler HPPTS 配置教程「建议收藏」 思路:先清除,之前的网站证书再安装fiddler的证书。具体看下面:https://www.cnblogs.com/xiao-qing/p/7771289.html转载于:https://www.cnblogs.com/limi2019/articles/11233846.html…

  • 活动目录的介绍:深入浅出Active Directory系列(一)

    活动目录的介绍:深入浅出Active Directory系列(一)

发表回复

您的电子邮箱地址不会被公开。

关注全栈程序员社区公众号