服务器之家

服务器之家 > 正文

Python 做曲线拟合和求积分的方法

时间:2021-05-09 00:43     来源/作者:vellerzheng

这是一个由加油站油罐传感器测量的油罐高度数据和出油体积,根据体积和高度的倒数,用截面积来描述油罐形状,求出拟合曲线,再用标准数据,求积分来验证拟合曲线效果和误差的一个小项目。 主要的就是首先要安装Anaconda  python库,然后来运用这些数学工具。

?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
###最小二乘法试验###
import numpy as np
import pymysql
from scipy.optimize import leastsq
from scipy import integrate
###绘图,看拟合效果###
import matplotlib.pyplot as plt
from sympy import *
 
 
path='E:\PythonProgram\oildata.txt'
 
lieh0 =[]   #初始第一列,油管高度
liev1 =[]   #初始第二列,油枪记录的体积
 
h_median =[]  # 高度相邻中位数
h_subtract =[]   #相邻高度作差
v_subtract =[]   #相邻体积作差
select_h_subtr =[]   #筛选的高度作差 ########
select_v_subtr =[]   #筛选的体积作差
 
VdtH=[]      #筛选的V 和 t 的 倒数。
 
def loadData(Spath,lie0,lie1):
 with open(Spath,'r') as f0:
   for i in f0:
    tmp=i.split()
    lie0.append(float(tmp[0]))
    lie1.append(float(tmp[2]))
 print ("source lie0",len(lie0))
 
 
def connectMySQL():
 db = pymysql.connect(host='10.**.**.**', user='root', passwd='***', db="zabbix", charset="utf8") # 校罐
 cur = db.cursor()
 
 try:
  # 查询
  cur.execute("SELECT * FROM station_snapshot limit 10 ")
  for row in cur.fetchall():
   # print(row)
   id = row[0]
   snapshot_id = row[1]
   DateTime = row[13]
   attr1V = row[5]
   attr2H = row[6]
   print("id=%d ,snapshot_id=%s,DateTime=%s,attr1V =%s, attr2H=%s ",
     (id, snapshot_id, DateTime, attr1V, attr2H))
 except:
  print("Error: unable to fecth data of station_stock")
 
 try:
  cur.execute("SELECT * FROM can_stock limit 5");
  for row in cur.fetchall():
   # print(row)
   stockid = row[0]
   stationid = row[1]
   DateTime = row[4]
   Volume = row[5]
   Height = row[8]
   print("stockid=%d ,stationid=%s,DateTime=%s,Volume =%f, Height=%f ",
     (stockid, stationid, DateTime, Volume, Height))
 except:
  print("Error: unable to fecth data of can_snapshot")
 
 cur.close()
 db.close()
 
 
def formatData(h_med,h_subtr,v_subtr):
 lh0 = lieh0[:]
 del lh0[0]
 print("lh0 size(): ",len(lh0))
 
 lh1 =lieh0[:]
 del lh1[len(lh1)-1]
 
 print("lh1 size() : ",len(lh1))
 
 lv0 =liev1[:]
 del lv0[0]
 #print (liev1)
 print ("Souce liev1 size() : ",len(liev1))
 print ("lv1 size() :",len(lv0))
 """
 lv1 =liev1[:]
 del lv1[len(lv1)-1]
 print("lv1 size(): ",len(lv1))
 """
 h_med[:] =[(x+y)/2 for x,y in zip(lh0,lh1)]  ###采样点(Xi,Yi)###
 print("h_med size() : ", len(h_med))
 
 h_subtr[:] = [(y-x) for x,y in zip(lh0,lh1)]
 print("h_subtr size() : ", len(h_subtr))
 # v_subtr[:] = [(y-x) for x,y in zip(lv0,lv1)]
 v_subtr[:] = lv0
 print("v_subtr size() : ", len(v_subtr))
 
 
def removeBadPoint(h_med,h_sub,v_sub):
 for val in h_sub:
  position=h_sub.index(val)
  if 0.01 > val > -0.01:
   del h_sub[position]
   del h_med[position]
   del v_sub[position]
 v_dt_h_ay = [(y/x) for x, y in zip(h_sub, v_sub)]
 return v_dt_h_ay
 
 
 
def selectRightPoint(h_med,h_subtr,v_dt_h_ay):
 for val in v_dt_h_ay:
  pos = v_dt_h_ay.index(val)
  if val > 20 :
   del v_dt_h_ay[pos]
   del h_med[pos]
   del h_subtr[pos]
 for val in v_dt_h_ay:
  ptr = v_dt_h_ay.index(val)
  if val < 14:
   del v_dt_h_ay[ptr]
   del h_med[ptr]
   del h_subtr[ptr]
 
 
def writeFile(h_mp, v_dt_h):
 s='\n'.join(str(num)[1:-1] for num in h_mp)
 v='\n'.join(str(vdt)[1:-1] for vdt in v_dt_h)
 open(r'h_2.txt','w').write(s)
 open(r'v_dt_h.txt','w').write(v)
 print("write h_median: ",len(h_mp))
 # print("V_dt also is (y-x) : ",v_dt_h,end="\n")
 print("Write V_dt_h : ",len(v_dt_h))
# file=open('data.txt','w')
# file.write(str(h_mp));
# file.close
 
 
def integralCalculate(coeff,xspace):
 vCalcute =[]
 x = Symbol('x')
 a, b, c, d = coeff[0]
 y = a * x ** 3 + b * x ** 2 + c * x + d
 i=0
 while (i< len(xspace)-1) :
  m = integrate(y, (x, xspace[i], xspace[i+1]))
  vCalcute.append(abs(m))
  i=i+1
 
 print("求导结果:",vCalcute)
 print("求导长度 len(VCalcute): ",len(vCalcute))
 return vCalcute
 
 ###需要拟合的函数func及误差error###
 
def func(p,x):
 a,b,c,d=p
 return a*x**3+b*x**2+c*x+d #指明待拟合的函数的形状,设定的拟合函数。
 
def error(p,x,y):
 return func(p,x)-y #x、y都是列表,故返回值也是个列表
 
def leastSquareFitting(h_mp,v_dt_hl):
 p0=[1,2,6,10#a,b,c 的假设初始值,随着迭代会越来越小
 #print(error(p0,h_mp,v_dt_h,"cishu")) #目标是让error 不断减小
 #s="Test the number of iteration" #试验最小二乘法函数leastsq得调用几次error函数才能找到使得均方误差之和最小的a~c
 Para=leastsq(error,p0,args=(h_mp,v_dt_hl)) #把error函数中除了p以外的参数打包到args中
 a,b,c,d=Para[0]   #leastsq 返回的第一个值是a,b,c 的求解结果,leastsq返回类型相当于c++ 中的tuple
 print(" a=",a," b=",b," c=",c," d=",d)
 plt.figure(figsize=(8,6))
 plt.scatter(h_mp,v_dt_hl,color="red",label="Sample Point",linewidth=3) #画样本点
 x=np.linspace(200,2200,1000)
 y=a*x**3+b*x**2+c*x+d
 
 integralCalculate(Para,h_subtract)
 plt.plot(x,y,color="orange",label="Fitting Curve",linewidth=2) #画拟合曲线
 # plt.plot(h_mp, v_dt_hl,color="blue", label='Origin Line',linewidth=1) #画连接线
 plt.legend()
 plt.show()
 
def freeParameterFitting(h_mp,v_dt_hl):
 z1 = np.polyfit(h_mp, v_dt_hl, 6) # 第一个拟合,自由度为6
  # 生成多项式对象
 p1 = np.poly1d(z1)
 print("Z1:")
 print(z1)
 print("P1:")
 print(p1)
 print("\n")
 x = np.linspace(400, 1700, 1000)
 plt.plot(h_mp, v_dt_hl, color="blue", label='Origin Line', linewidth=1) # 画连接线
 plt.plot(x, p1(x), 'gv--', color="black", label='Poly Fitting Line(deg=6)', linewidth=1)
 plt.legend()
 plt.show()
 
def main():
 loadData(path, lieh0, liev1)
 connectMySQL() # 读取oildata数据库
 
 formatData(h_median, h_subtract, v_subtract)
 
 # 去除被除数为0对应的点,并得到v 和 h 求导 值的列表
 VdtH[:] = removeBadPoint(h_median, h_subtract, v_subtract)
 print("h_median1:", len(h_median))
 
 print("VdtH1 : ", len(VdtH))
 
 # 赛选数据,去除异常点
 selectRightPoint(h_median, h_subtract, VdtH)
 print("h_median2:", len(h_median))
 print("h_subtract: ", len(h_subtract))
 print("VdtH2 : ", len(VdtH))
 h_mp = np.array(h_median)
 v_dt_h = np.array(VdtH)
 
 writeFile(h_mp, v_dt_h)
 # 最小二乘法作图
 leastSquareFitting(h_mp, v_dt_h)
 # 多项式自由参数法作图
 freeParameterFitting(h_mp, v_dt_h)
 
if __name__ == '__main__':
 main()

以上这篇Python 做曲线拟合和求积分的方法就是小编分享给大家的全部内容了,希望能给大家一个参考,也希望大家多多支持服务器之家。

原文链接:https://blog.csdn.net/vellerzheng/article/details/71544511

相关文章

热门资讯

2020微信伤感网名听哭了 让对方看到心疼的伤感网名大全
2020微信伤感网名听哭了 让对方看到心疼的伤感网名大全 2019-12-26
yue是什么意思 网络流行语yue了是什么梗
yue是什么意思 网络流行语yue了是什么梗 2020-10-11
背刺什么意思 网络词语背刺是什么梗
背刺什么意思 网络词语背刺是什么梗 2020-05-22
苹果12mini价格表官网报价 iPhone12mini全版本价格汇总
苹果12mini价格表官网报价 iPhone12mini全版本价格汇总 2020-11-13
2021德云社封箱演出完整版 2021年德云社封箱演出在线看
2021德云社封箱演出完整版 2021年德云社封箱演出在线看 2021-03-15
返回顶部