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

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

###最小二乘法试验###
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()

来源：https://blog.csdn.net/vellerzheng/article/details/71544511