你希望通过几种常见算法的实现,了解python在数学建模中的能力。
python除了丰富的原生数据结构外,拥有强大的第三方软件包支持,例如矩阵运算库Numpy,数据处理库Pandas、机器学习库Sklearn、深度学习库Tenserflow&Pytorch、科学计算库Scipy、图形绘制库matplotlib、网络算法库Networkx。此外几乎针对任何领域,都有第三方软件包的支持,这归功于python优秀的社区。使用者需要使用好pip这一软件包管理工具,发掘前人造好的轮子,尽量减少自己编程的难度。我们将在后面的问题讨论中介绍以下几种常用数学建模算法的python实现:
4.单源多宿最短路算法
我们的重点在于代码实现而非数学推导
1.数据拟合算法
我们这里介绍通过最小二乘法拟合线性函数
-
#我们使用最小二乘法拟合一个三次函数,选取了5个参数 -
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()

-
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()

-
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])})#定义函数的多个约束条件 -
res=minimize(func,[-1.0,1.0],constraints=cons,method="SLSQP",options={"disp":True}) -
print(res)

-
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