Python日记（2）——反应速率的计算

Python日记（2）——反应速率的计算

每天做一个Python小练习，顺便记录一些小技巧。
化工数值计算实验的实验二：1，3-丁二烯气相二聚反应的速率计算
一、目的和要求
1.掌握常用的插值方法
2.掌握数值微分的计算方法
二、问题提出
1.函数插值
函数是工程实践中重要的工具，它能够描述变量之间的关系。但是实际应用
当中，我们很难使用解析表达式描述两变量之间的关系。我们能够得到的，只是
从实验测量得到的离散的数据，即函数表。当我们需要使用的数据点不在函数表
中时，我们就需要进行插值。
2.数值微分
有时候，需要从实验得到的有限的数据中推断出变量之间的相对变化率。但
是，仅从有限个数据，不可能通过解析的方法精确求解出其导数，只能利用数值
方法计算。数值微分与插值和拟合是密不可分的。在进行数值计算时，往往需要
插值或拟合来增加数据点的密度，提高精度。
3.二聚反应速率的计算

t/min	p/mmHg
0	632.0
5	611.0
10	592.0
15	573.5
20	558.5
25	545.0
30	533.0
35	523.0
40	514.0
45	505.0
50	497.0
55	490.0
60	484.0
65	478.5
70	473.0
75	468.0
80	463.0
90	458.0
95	453.0

丁二烯可发生二聚反应：
2C4H6 → (C4H6)2
若初始压力为p0，在t时刻丁二烯的分压为p1，二聚体的分压为p2，则它们存
在关系p1+ 2p2= p0。该反应的速率可表示为
-1/2dp1/dt = dp2/dt
也就是“虽计+算反应速率，需要知道P或p2与时间的函数关系，但是实验只能测量总压，而且实验测量只能得到有限的离散的数据，因此需要从有限的数据中计算分压对时间的导数，从而得到反应速率。
三、实验内容
1.换算p的单位
根据 1kPa = 7.5mmHg

def p_convert(p):
    for i in p:
        index_num = p.index(i)
        i /= 7.5
        p[index_num] = round(i, 4)
    return p

2.因为二聚体的分压等于Δp，所以利用循环求Δp并返回便是p1的分压

def p_delta(p):
    p1 = list()
    for i in range(0, len(p)):
        if i+1 == len(p):
            break
        else:
            pl = round(p[0] - p[i])
            p1.append(pl)
    return p1

3.曲线插值
（1）设置x的点集
linspace(begin, end, num)
begin为左端点，end为右端点，num为元素个数

def generate_data(begin, end, num):
    """ 产生x点集 """
    x = np.linspace(begin, end, num)
    return x

（1）线性插值
interp1d(x, y, kind)
x, y为已知数据点，kind为可选类型，可以是字符串或整数

def interpolate_linear(x, y):
    """ 线性插值 """
    f_linear = interpolate.interp1d(x, y)
    return f_linear

(3)求得样条曲线或导数，默认为0阶导数。
scipy.interpolate.splrep(x, y, w=None, xb=None, xe=None, k=3, task=0, s=None, t=None, full_output=0, per=0, quiet=1)
查找一维曲线的B样条曲线表示。给定数据点集，确定区间上度k的平滑样条近似。
（各个参数的表达含义我还没学会）
x, y为已知数据点，x_new就是新的细分后的x点集，元素个数要与差值后的个数保持一直。

def interpolate_b_spline(x, y, x_new, der=0):
    """ B样条曲线插值或者导数， 默认der=0 """
    tck = interpolate.splrep(x, y)
    y_bspline = interpolate.splev(x_new, tck, der=der)
    return y_bspline

(4)主程序

def num_interpolate(begin, end, p, p0):
    pt_x1 = generate_data(begin=begin, end=end, num=len(p))
    pt_x2 = generate_data(begin=begin, end=end, num=len(p0))
    pt_y1 = p
    pt_y2 = p0

    interpolate_x = generate_data(begin=begin, end=end, num=900)

    f_linear1 = interpolate_linear(pt_x1, pt_y1)
    f_linear2 = interpolate_linear(pt_x2, pt_y2)

    y_bspline1 = interpolate_b_spline(pt_x1, pt_y1, interpolate_x, der=0)
    y_bspline1_derivative = interpolate_b_spline(pt_x1, pt_y1, interpolate_x, der=1)
    y_bspline2 = interpolate_b_spline(pt_x2, pt_y2, interpolate_x, der=0)
    y_bspline2_derivative = interpolate_b_spline(pt_x1, pt_y1, interpolate_x, der=1)

    fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(15, 20))
    ax1.plot(pt_x1, pt_y1, 'o', label='Origin Data')        # 原始数据点
    ax1.plot(interpolate_x, y_bspline1, color='blue')       # 插值后的平滑曲线
    ax1.plot(interpolate_x, f_linear1(interpolate_x), 'b--', label=u'丁二烯的分压')
    ax1.plot(interpolate_x, y_bspline2, color='red')
    ax1.plot(interpolate_x, f_linear2(interpolate_x), 'r--', label=u'二聚体的分压')
    ax1.legend()

    ax2.plot(interpolate_x, -y_bspline1_derivative, label='速率与时间的关系')
    ax3.plot(interpolate_x, -y_bspline2_derivative, label='速率与时间的关系')
    plt.show()


if __name__ == '__main__':
    p = [632.0, 611.0, 592.0, 573.5, 558.5, 545.0, 533.0, 523.0, 514.0, 505.0, 497.0, 490.0, 484.0, 478.5, 473.0, 468.0,
         463.0, 458.0, 453.0]

    num_interpolate(
        begin=0,
        end=90,
        p=p_convert(p),
        p0=p_delta(p)
    )

(4)结果图运行代码得

总代码如下：

import numpy as np
from matplotlib import pyplot as plt
from scipy import interpolate
import matplotlib


matplotlib.rcParams['font.sans-serif'] = ['Microsoft YaHei']


def p_convert(p):
    for i in p:
        index_num = p.index(i)
        i /= 7.5
        p[index_num] = round(i, 4)
    return p


def p_delta(p):
    p1 = list()
    for i in range(0, len(p)):
        if i+1 == len(p):
            break
        else:
            pl = round(p[0] - p[i])
            p1.append(pl)
    return p1


def generate_data(begin, end, num):
    """ 产生x点集 """
    x = np.linspace(begin, end, num)
    return x


def interpolate_linear(x, y):
    """ 线性插值 """
    f_linear = interpolate.interp1d(x, y)
    return f_linear


def interpolate_b_spline(x, y, x_new, der=0):
    """ B样条曲线插值或者导数， 默认der=0 """
    tck = interpolate.splrep(x, y)
    y_bspline = interpolate.splev(x_new, tck, der=der)
    return y_bspline


def num_interpolate(begin, end, p, p0):
    pt_x1 = generate_data(begin=begin, end=end, num=len(p))
    pt_x2 = generate_data(begin=begin, end=end, num=len(p0))
    pt_y1 = p
    pt_y2 = p0

    interpolate_x = generate_data(begin=begin, end=end, num=900)

    f_linear1 = interpolate_linear(pt_x1, pt_y1)
    f_linear2 = interpolate_linear(pt_x2, pt_y2)

    y_bspline1 = interpolate_b_spline(pt_x1, pt_y1, interpolate_x, der=0)
    y_bspline1_derivative = interpolate_b_spline(pt_x1, pt_y1, interpolate_x, der=1)
    y_bspline2 = interpolate_b_spline(pt_x2, pt_y2, interpolate_x, der=0)
    y_bspline2_derivative = interpolate_b_spline(pt_x1, pt_y1, interpolate_x, der=1)

    fig, (ax1, ax2, ax3) = plt.subplots(3, 1)
    ax1.plot(pt_x1, pt_y1, 'o', label='Origin Data')        # 原始数据点
    ax1.plot(interpolate_x, y_bspline1, color='blue')       # 插值后的平滑曲线
    ax1.plot(interpolate_x, f_linear1(interpolate_x), 'b--', label=u'丁二烯的分压')
    ax1.plot(interpolate_x, y_bspline2, color='red')
    ax1.plot(interpolate_x, f_linear2(interpolate_x), 'r--', label=u'二聚体的分压')
    ax1.legend()

    ax2.plot(interpolate_x, -y_bspline1_derivative, label='速率与时间的关系')
    ax3.plot(interpolate_x, -y_bspline2_derivative, label='速率与时间的关系')
    plt.show()


if __name__ == '__main__':
    p = [632.0, 611.0, 592.0, 573.5, 558.5, 545.0, 533.0, 523.0, 514.0, 505.0, 497.0, 490.0, 484.0, 478.5, 473.0, 468.0,
         463.0, 458.0, 453.0]

    num_interpolate(
        begin=0,
        end=90,
        p=p_convert(p),
        p0=p_delta(p)
    )

原方法转载自：原文链接：https://blog.csdn.net/comedate/article/details/109532850