如何使用python完成數學建模常見算法!

在數學建模中主流的編程語言是MATLAB,但隨著python/R中數學軟件包的不斷完善,熟悉這兩種編程語言的同學也可以快速數學建模的編程環節。後麵我們將介紹幾種常見數學建模算法的python實現,旨在展示python在本領域的強大威力。

1、問題描述

你希望通過幾種常見算法的實現,了解python在數學建模中的能力。

2、解決方案

python除了豐(feng) 富的原生數據結構外,擁有強大的第三方軟件包支持,例如矩陣運算庫Numpy,數據處理庫Pandas、機器學習(xi) 庫Sklearn、深度學習(xi) 庫Tenserflow&Pytorch、科學計算庫Scipy、圖形繪製庫matplotlib、網絡算法庫Networkx。此外幾乎針對任何領域,都有第三方軟件包的支持,這歸功於(yu) python優(you) 秀的社區。使用者需要使用好pip這一軟件包管理工具,發掘前人造好的輪子,盡量減少自己編程的難度。我們(men) 將在後麵的問題討論中介紹以下幾種常用數學建模算法的python實現:

1.數據擬合算法
2.插值算法
3.線性規劃算法

4.單源多宿最短路算法

3、問題討論

我們(men) 的重點在於(yu) 代碼實現而非數學推導

1.數據擬合算法

我們(men) 這裏介紹通過最小二乘法擬合線性函數

  1. #我們(men) 使用最小二乘法擬合一個(ge) 三次函數,選取了5個(ge) 參數

  2. import numpy as np

  3. import matplotlib.pyplot as plt

  4. SAMPLE_NUM = 100

  5. M = 5

  6. x = np.arange(0, SAMPLE_NUM).reshape(SAMPLE_NUM, 1) / (SAMPLE_NUM - 1) * 10

  7. y = 2*x3+3*x2+x+1

  8. plt.plot(x, y, 'bo')

  9. X = x

  10. for i in range(2, M + 1):

  11.    X = np.column_stack((X, pow(x, i)))

  12. X = np.insert(X, 0, [1], 1)

  13. W=np.linalg.inv((X.T.dot(X))).dot(X.T).dot(y)

  14. y_estimate = X.dot(W)

  15. plt.plot(x, y_estimate, 'r')

  16. plt.show()

國賽必備 | 了解如何使用python完成數學建模常見算法!

2.插值算法
我們使用幾種常見的插值函數擬合上例中的三次函數
  1. import numpy as np

  2. from scipy import interpolate

  3. import pylab as pl

  4. x=np.linspace(0,10,11)

  5. y=2*x3+3*x2+x+1

  6. xInset=np.linspace(0,10,101)

  7. pl.plot(x,y,"ro")

  8. for kind in["nearest","zero","slinear","quadratic","cubic"]:

  9.    f=interpolate.interp1d(x,y,kind=kind)

  10.    y_estimate=f(xInset)

  11.    pl.plot(xInset,y_estimate,label=str(kind))

  12. pl.legend(loc="lower right")

  13. pl.show()

國賽必備 | 了解如何使用python完成數學建模常見算法!

3.線性規劃算法
我們可以使用scipy庫對目標函數進行線性規劃
  1. import numpy as np

  2. from scipy.optimize import minimize

  3. def func(x):

  4. return(2*x[0]*x[1]+2*x[0]-x[0]2+2*x[1]2+np.sin(x[0]))

  5. cons=({"type":"eq","fun":lambda x:np.array([x[0]3-x[1]]),"jac":lambda x:np.array([3*(x[0]2),-1.0])},{"type":"ineq","fun":lambda x:np.array([x[1]-1]),"jac":lambda x:np.array([0,1])})#定義(yi) 函數的多個(ge) 約束條件

  6. res=minimize(func,[-1.0,1.0],constraints=cons,method="SLSQP",options={"disp":True})

  7. print(res)

注意這裏求解的為目標函數最小值,如果我們需要求解最大值則將func取負數即可。輸出內容如圖

國賽必備 | 了解如何使用python完成數學建模常見算法!

4.單源多宿最短路算法
我們介紹一下基於堆優化的dijkstra算法,這裏的堆可以使用python內置的PriorityQueue實現。我們這裏給出一個簡單的堆實現:
  1. classDisNode:

  2. def __init__(self,node,dis):

  3.        self.node=node

  4.        self.dis=dis

  5. def __lt__(self, other):

  6. return self.dis<other.dis

  7. classDisPath:

  8. def __init__(self,end):

  9.        self.end=end

  10.        self.path=[self.end]

  11.        self.dis=0

  12. def __str__(self):

  13.        nodes=self.path.copy()

  14. return"->".join(list(map(str,nodes)))+" "+str(self.dis)

  15. classHeap:

  16. def __init__(self):

  17.        self.size=0

  18.        self.maxsize=10000

  19.        self.elements=[0]*(self.maxsize+1)

  20. def isEmpty(self):

  21. return self.size==0

  22. def insert(self,value):

  23. if self.isEmpty():

  24.            self.elements[1]=value

  25. else:

  26.            index=self.size+1

  27. while(index!=1and value<self.elements[index//2]):

  28.                self.elements[index]=self.elements[index//2]

  29.                index=index//2

  30.            self.elements[index]=value

  31.        self.size+=1

  32. def pop(self):

  33.        deleteElement=self.elements[1]

  34.        self.elements[1]=self.elements[self.size]

  35.        self.size-=1

  36.        temp=self.elements[1]

  37.        parent,child=1,2

  38. while(child<=self.size):

  39. if child<self.size and self.elements[child]>self.elements[child+1]:

  40.                child+=1

  41. if temp<self.elements[child]:

  42. break

  43. else:

  44.                self.elements[parent]=self.elements[child]

  45.                parent=child

  46.                child*=2

  47.        self.elements[parent]=temp

  48. return deleteElement

  49. defDijkstraWithHeap(nodes,start,GetNeighbors):

  50.    dis=defaultdict(int)

  51.    paths=defaultdict(DisPath)

  52.    heap=Heap()

  53.    visit=set()

  54. for node in nodes:

  55.        dis[node]=sys.maxsize

  56.        paths[node]=DisPath(node)

  57.    dis[start]=0

  58.    heap.insert(DisNode(start,0))

  59. while(not heap.isEmpty()):

  60.        now=heap.pop().node

  61. if now in visit:

  62. continue

  63.        visit.add(now)

  64.        paths[now].dis=dis[now]

  65. for edge inGetNeighbors(now):

  66.            end=edge.End

  67. if dis[now]+edge.value<dis[end]:

  68.                dis[end]=dis[now]+edge.value

  69.                paths[end].path=paths[now].path+[end]

  70.                heap.insert(DisNode(end,dis[end]))

  71. return paths

通過堆優化的dijkstra算法的時間複雜度最低可以達到O(nlogm),上文給出的代碼的時間複雜度為O(mlogm),但是勝在編碼簡單。此外還可以使用Networkx庫進行求解,這裏不再介紹。

【競賽報名/項目谘詢+微信:mollywei007】

上一篇

美國中學生7大寫作競賽及寫作活動盤點

下一篇

2022 iGEM國際基因工程機器競賽組隊報名中

你也可能喜歡

  • 暫無相關文章!

評論已經被關(guan) 閉。

插入圖片
返回頂部