Python 做曲线拟合和求积分

来源：互联网发布：图像相似度算法编辑：程序博客网时间：2024/05/29 14:55
这是一个由加油站油罐传感器测量的油罐高度数据和出油体积，根据体积和高度的倒数，用截面积来描述油罐形状，求出拟合曲线，再用标准数据，求积分来验证拟合曲线效果和误差的一个小项目。主要的就是首先要安装Anaconda python库，然后来运用这些数学工具。
###最小二乘法试验###import numpy as npimport pymysqlfrom scipy.optimize import leastsqfrom scipy import integrate###绘图，看拟合效果###import matplotlib.pyplot as pltfrom 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_aydef 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.closedef 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()
0 0