數學建模比賽必備：用python完成數學建模常見算法

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

1、問題描述

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

2、解決(jue) 方案

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) 這裏介紹通過最小二乘法擬合線性函數

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

import numpy as np

import matplotlib.pyplot as plt

SAMPLE_NUM = 100

M = 5

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

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

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

X = x

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

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

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

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

y_estimate = X.dot(W)

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

plt.show()

2.插值算法

我們(men) 使用幾種常見的插值函數擬合上例中的三次函數

import numpy as np

from scipy import interpolate

import pylab as pl

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

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

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

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

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

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

y_estimate=f(xInset)

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

pl.legend(loc="lower right")

pl.show()

3.線性規劃算法

我們(men) 可以使用scipy庫對目標函數進行線性規劃

import numpy as np

from scipy.optimize import minimize

def func(x):

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

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) 約束條件

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

print(res)

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

4.單源多宿最短路算法

我們(men) 介紹以下基於(yu) 堆優(you) 化的dijkstra算法，這裏的堆可以使用python內(nei) 置的PriorityQueue實現。我們(men) 這裏給出一個(ge) 簡單的堆實現：

classDisNode:

def __init__(self,node,dis):

self.node=node

self.dis=dis

def __lt__(self, other):

return self.dis<other.dis

classDisPath:

def __init__(self,end):

self.end=end

self.path=[self.end]

self.dis=0

def __str__(self):

nodes=self.path.copy()

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

classHeap:

def __init__(self):

self.size=0

self.maxsize=10000

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

def isEmpty(self):

return self.size==0

def insert(self,value):

if self.isEmpty():

self.elements[1]=value

else:

index=self.size+1

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

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

index=index//2

self.elements[index]=value

self.size+=1

def pop(self):

deleteElement=self.elements[1]

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

self.size-=1

temp=self.elements[1]

parent,child=1,2

while(child<=self.size):

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

child+=1

if temp<self.elements[child]:

break

else:

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

parent=child

child*=2

self.elements[parent]=temp

return deleteElement

defDijkstraWithHeap(nodes,start,GetNeighbors):

dis=defaultdict(int)

paths=defaultdict(DisPath)

heap=Heap()

visit=set()

for node in nodes:

dis[node]=sys.maxsize

paths[node]=DisPath(node)

dis[start]=0

heap.insert(DisNode(start,0))

while(not heap.isEmpty()):

now=heap.pop().node

if now in visit:

continue

visit.add(now)

paths[now].dis=dis[now]

for edge inGetNeighbors(now):

end=edge.End

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

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

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

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

return paths

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