2021高教社杯全國大學生數學建模

B 題 乙醇偶合制備 C4 烯烴

C4 烯烴廣泛應用于化工產品及醫藥的生產，乙醇是生產制備 C4 烯烴的原料，在制備程序中，催化劑組合（即：Co 負載量、Co/SiO2 和 HAP 裝料比、乙醇濃度的組合）與溫度對 C4 烯烴的選擇性和 C4 烯烴收率將產生影響（名詞解釋見附錄），
因此通過對催化劑組合設計，探索乙醇催化偶合制備 C4 烯烴的工藝條件具有非常重要的意義和價值，
某化工實驗室針對不同催化劑在不同溫度下做了一系列實驗，結果如附件 1 和附件 2 所示，請通過數學建模完成下列問題：
(1) 對附件 1 中每種催化劑組合，分別研究乙醇轉化率、C4 烯烴的選擇性與溫度的關系，并對附件 2 中 350 度時給定的催化劑組合在一次實驗不同時間的測驗結果進行分析，
(2) 探討不同催化劑組合及溫度對乙醇轉化率以及 C4 烯烴選擇性大小的影響，
(3) 如何選擇催化劑組合與溫度，使得在相同實驗條件下 C4 烯烴收率盡可能高，若使溫度低于 350 度，又如何選擇催化劑組合與溫度，使得 C4 烯烴收率盡可能高，
(4) 如果允許再增加 5 次實驗，應如何設計，并給出詳細理由，

問題一程式代碼：

from mpl_toolkits import mplot3d
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
m=200; n=300
x=np.linspace(-6, 6, m); y=np.linspace(-8, 8, n);
x2, y2 = np.meshgrid(x, y)
x3=np.reshape(x2,(1,-1)); y3=np.reshape(y2, (1,-1))
xy=np.vstack((x3,y3))
def Pfun(t, m1, m2, s):
    return np.exp(-((t[0]-m1)**2+(t[1]-m2)**2)/(2*s**2))
z=Pfun(xy, 1, 2, 3); zr=z+0.2*np.random.normal(size=z.shape) #噪聲資料
popt, pcov=curve_fit(Pfun, xy, zr)   #擬合引數
print("三個引數的擬合值分別為：",popt)
zn=Pfun(xy, *popt)  #計算擬合函式的值
zn2=np.reshape(zn, x2.shape)
plt.rc('font',size=16)
ax=plt.axes(projection='3d') #創建一個三維坐標軸物件
ax.plot_surface(x2, y2, zn2,cmap='gist_rainbow')
plt.savefig("figure7_10.png", dpi=500); plt.show()

import matplotlib.pyplot as plt
 
"""
@brief:   計算n階差商 f[x0, x1, x2 ... xn]
@param:   xi   所有插值節點的橫坐標集合                                                        o
@param:   fi   所有插值節點的縱坐標集合                                                      /   \
@return:  回傳xi的i階差商(i為xi長度減1)                                                     o     o
@notice:  a. 必須確保xi與fi長度相等                                                        / \   / \
          b. 由于用到了遞回，所以留意不要爆堆疊了.                                           o   o o   o
          c. 遞回減遞回(每層遞回包含兩個遞回函式), 每層遞回次數呈二次冪增長，總次數是一個滿二叉樹的所有節點數量(所以極易堆疊溢位)                                                                                    
"""
def diff_quo(xi = [], fi = []):
    if len(xi) > 2 and len(fi) > 2:
        return (diff_quo(xi[:len(xi) - 1], fi[:len(fi) - 1]) - diff_quo(xi[1:len(xi)], fi[1:len(fi)])) / float(xi[0] - xi[-1])
    return (fi[0] - fi[1]) / float(xi[0] - xi[1])
     
 
 
 
"""
@brief:  獲得Wi(x)函式;
         Wi的含義舉例 W1 = (x - x0); W2 = (x - x0)(x - x1); W3 = (x - x0)(x - x1)(x - x2)
@param:  i  i階(i次多項式)
@param:  xi  所有插值節點的橫坐標集合
@return: 回傳Wi(x)函式
"""
def get_Wi(i = 0, xi = []):
    def Wi(x):
        result = 1.0
        for each in range(i):
            result *= (x - xi[each])
        return result
    return Wi
     
     
     
     
"""
@brief: 獲得牛頓插值函式
@
"""
def get_Newton_inter(xi = [], fi = []):
    def Newton_inter(x):
        result = fi[0]
        for i in range(2, len(xi)):
            result += (diff_quo(xi[:i], fi[:i]) * get_Wi(i-1, xi)(x))
        return result
    return Newton_inter
     
         
 
"""
demo:
"""
if __name__ == '__main__':  
 
    ''' 插值節點, 這里用二次函式生成插值節點，每兩個節點x軸距離位10 '''
    sr_x = [i for i in range(-50, 51, 10)]
    sr_fx = [i**2 for i in sr_x] 
     
    Nx = get_Newton_inter(sr_x, sr_fx)            # 獲得插值函式
     
    tmp_x = [i for i in range(-50, 51)]          # 測驗用例
    tmp_y = [Nx(i) for i in tmp_x]               # 根據插值函式獲得測驗用例的縱坐標
        
    ''' 畫圖 '''
    plt.figure("aaa")
    ax1 = plt.subplot(111)
    plt.sca(ax1)
    plt.plot(sr_x, sr_fx, linestyle = '', marker='o', color='b')
    plt.plot(tmp_x, tmp_y, linestyle = '--', color='r')
    plt.show()

解題關鍵

print("聯系qq：3506426881")