这学期选了数学建模课，因为上学期学了MATLAB，这学期尝试使用Python完成数学建模。
Python的基本语法其实很简单，这里推荐菜鸟教程https://www.runoob.com/python/python-basic-syntax.html和廖雪峰的python教程https://www.liaoxuefeng.com/wiki/1016959663602400
python的强大在于他的各种包，难也难在各种包。要想熟练使用各种包中的各种函数还是有一定难度的，有时候不知道为什么就掉坑里了。
昨天用python写了如下几个问题：
1. 建立M-文件: 已知函数 
计算f (-1), f (0.5), f (1.5)，并作出该函数的曲线图。
2. 编写利用顺序Guass 消去法求方程组解的M-函数文件,并计算方程组
的解
3. 编写“商人们安全过河”的Matlab 程序
4. 编写“人口预报问题”的Matlab 程序
上面的matlab都换成python
1. 第一个问题比较基础，写的时候只遇到了一个小问题。
1 import numpy as np 2 import matplotlib.pyplot as plt 3 def myfunc(x): 4 if 0 > x >= -1: 5 return x+1 6 elif 1 > x >= 0: 7 return 1 8 elif 2 > x >= 1: 9 return x**2 10 11 x = np.linspace(-1, 2, 300) 12 y = np.array([myfunc(t) for t in x]) 13 14 y1 = myfunc(-1) 15 y2 = myfunc(0.5) 16 y3 = myfunc(1.5) 17 print("f(-1) = %.2f f(0.5) = %.2f f(1.5) = %.2f" %(y1, y2, y3)) 18 plt.figure() 19 plt.plot(x, y) 20 plt.show()
print函数的格式化输出！我看教程的时候觉得和C语言一样的就直接略过了，结果自己写的时候才发现python的格式化输出是有一些不同的，主要就是后面变量的表示用%和前面的字符串分割而不是C语言中的 ，
另一种格式化字符串的方法是使用字符串的format()方法，它会用传入的参数依次替换字符串内的占位符{0}、{1}……，不过这种方式写起来比%要麻烦得多：
1 ‘Hello, {0}, 成绩提升了 {1:.1f}%‘.format(‘小明‘, 17.125) 2 ‘Hello, 小明, 成绩提升了 17.1%‘
2. 高斯消元法
1 import matplotlib.pyplot as plt 2 import numpy as np 3 from pylab import mpl 4 import math 5 # step0 消元 6 def step0(matrix): 7 row = matrix.shape[0] 8 9 # 保证主元为一 或者主元所在行全为 0 10 for i in range(0, row): 11 if not matrix[i, i]: 12 for j in range(i + 1, row): 13 if matrix[j, i]: 14 matrix[[i, j], :] = matrix[[j,i], :] 15 break 16 17 # 开始消元 18 for i in range(0, row - 1): # 以这些行的主元作为参照依次消除主元以下元素 19 for j in range(i + 1, row): 20 matrix[j, :] = matrix[j, :] - matrix[i, :]/matrix[i, i]*matrix[j,i] 21 22 return matrix 23 24 # step1 回代 25 def step1(matrix): 26 row = matrix.shape[0] 27 # 从倒数第二行开始消元 28 for i in range(row - 2, -1, -1): 29 for j in range(i + 1, row): 30 matrix[i, :] = matrix[i, :] - matrix[j, :]/matrix[j, j]*matrix[i,j] 31 32 return matrix 33 34 # 高斯消元法 35 def Gauss(matrix): 36 37 new = step1(step0(matrix)) 38 row = new.shape[0] 39 col = new.shape[1] 40 x = [0 for k in range (col-1) ] 41 for i in range(0,row): 42 43 x[i] = new[i,col-1]/new[i,i] 44 45 46 return x 47 48 paremeters=np.matrix(‘1,1,-1,1; 1,2, -2, 0; -2,1,1,1‘) 49 a = Gauss(paremeters) 50 print(a)
这里借鉴了一下网上的代码。主函数是全部自己写的。遇到的问题是：for循环的范围问题。
这里要强调两点
1）range（n，m）表示的是n到m-1（包括）的数，默认步长为1，可以在n，m中间加入参数自定义步长。
2） np.matrix定义的矩阵，行号和列号是从0开始的；python自带的list，索引也是从0开始的。
把握好这两点就可以避免掉坑里
3. 商人过河问题
1 #允许状态集合，例num=3 2 #S={(x,y)|x=0,y=0,1,2,3;x=3,y=0,1,2,3;x=y=1,2} x是此岸的商人数，y是此岸的仆人数 3 #允许决策集合，例boat_limit=2 4 #D={(u,v)|1<=u+v<=2,u,v=0,1,2} u是撘载的商人数，v是搭载的仆人数 5 def cross_river(person, boat): 6 num = person 7 boat_limit= boat 8 temp=[]#允许状态集合 9 for i in range(0,num+1): 10 if i==0 or i==num: 11 for j in range(0,num+1): 12 temp.append((i,j)) 13 else: 14 temp.append((i,i)) 15 S=set(temp) 16 D=[]#允许决策集合 17 for u in range(0,boat_limit+1): 18 for v in range(0,boat_limit+1): 19 if u+v>=1 and u+v<=boat_limit : 20 D.append((u,v)) 21 start=(num,num)#起始 22 end=(0,0)#目标 23 queue=[] 24 queue.append((0,start)) #前面的元素如果是0，说明是船在此岸，是1，说明船在对岸 25 step_dict={} 26 flag=0 27 finish=[] 28 while len(queue)!=0: 29 q_pop=queue.pop(0) 30 if q_pop[0]==0: 31 for x in D: 32 temp_s=(q_pop[1][0]-x[0],q_pop[1][1]-x[1]) 33 if temp_s not in S: 34 continue 35 if (1,temp_s) in step_dict: 36 continue 37 queue.append((1,temp_s)) 38 step_dict[(1,temp_s)]=q_pop 39 if temp_s==end: 40 flag=1 41 finish=(1,temp_s) 42 break 43 else: 44 for x in D: 45 temp_s=(q_pop[1][0]+x[0],q_pop[1][1]+x[1]) 46 if temp_s not in S: 47 continue 48 if (0,temp_s) in step_dict: 49 continue 50 queue.append((0,temp_s)) 51 step_dict[(0,temp_s)]=q_pop 52 if flag==1: 53 break 54 if flag==1: 55 print(‘该问题有解！最短路径：‘) 56 path=[] 57 path.append(finish) 58 while path[-1]!=(0,start): 59 path.append(step_dict[path[-1]]) 60 path.reverse() 61 real_path=list(map(str,path)) 62 for i in range(len(real_path)): 63 if i!=len(real_path)-1: 64 print(real_path[i] + ‘->‘) 65 else: 66 print(real_path[i]) 67 else: 68 print(‘该问题无解‘) 69 70 return None 71 72 73 cross_river(3,2)
这个问题简单来说： 三名商人各自带领一个随从渡河，一只小船只能容纳两人，在任意岸边，随从人数不能大于商人人数，不然就会杀人越货。
解决方法：用状态（变量）表示某一岸的人员状况；决策（变量）表示船上的人员状况，可以找出状态随决策变化的规律，在状态的允许变化范围内（即安全渡河条件），确定每一步的决策，达到渡河目标
4. 人口预测
（1） 假设一：指数增长模型，人口（相对）增长率r是常数
（2） 假设二：阻滞增长模型-logestic模型，自然资源和环境因素对人口的增长期阻滞作用，人口规模增大时，人口增长率降低,且假定r与人口x的关系是线性的；
也就是有两种公式来预测，首先经过数学运算得到两个函数表达式，用python分别定义函数，再用cruve_fit（）拟合出函数中未定的参数。详细的问题介绍和推导可以百度到。
1 import matplotlib.pyplot as plt 2 from pylab import mpl 3 import math 4 import numpy as np 5 import matplotlib.pyplot as plt 6 from scipy.optimize import curve_fit 7 8 t = [] 9 n = 0 10 for i in range(1790, 2010, 10): 11 12 t.append (i) 13 n = n+1 14 15 y = [3.9, 5.3, 7.2, 9.6, 12.9, 17.1, 23.2, 31.4, 16 38.6, 50.2, 62.9, 76.0,92.0,106.5, 123.2,131.7, 17 150.7, 179.3,204.0, 226.5, 251.4, 281.4] 18 19 20 #指数增长 21 def func(t,r): 22 return y[0]*np.exp(r*t) 23 24 t2 = np.array(t)-1790 25 y2 = np.array(y) 26 popt, pcov = curve_fit(func, t2, y2, p0 = [0.02]) 27 #如果直接这样拟合要给个初值p0，不然会拟合次数达到上限，拟合失败 28 print(popt) 29 r = popt[0] 30 31 #阻滞增长 32 def func2(t, r2, x_m): 33 return x_m/(1+ (((x_m/y[0])-1)*np.exp((-r2)*t))) 34 popt2, pcov2 = curve_fit(func2, t2, y2, p0 = [r, 400]) 35 #如果直接这样拟合要给个初值p0，不然会拟合结果和实际相差太多 36 print(popt2) 37 r2 = popt2[0] 38 x_m = popt2[1] 39 40 #画图 41 yvals = func(t2, r) 42 yvals2 = func2(t2,r2, x_m) 43 plt.figure() 44 plot1 = plt.plot(t2, y2, ‘s‘,label=‘original values‘) 45 plot2 = plt.plot(t2, yvals, ‘r‘,label=‘exp‘) 46 plot3 = plt.plot(t2, yvals2, ‘g‘,label=‘zuzhi‘) 47 48 plt.xlabel(‘t‘) 49 plt.ylabel(‘y‘) 50 plt.legend(loc=4) #指定legend的位置右下角 51 plt.title(‘curve_fit‘) 52 plt.show() 53 54 #预测2010年的人口 55 pred = func2(2010-1790, r2, x_m ) 56 print(pred)
这里两个函数都是直接把初始形式带进去拟合的，这会很难拟合，有时候会拟合次数达到上限，拟合失败，有时候拟合结果和实际相差太多；如果要直接用初始形式拟合，一定要给cruve_fit()传入p0参数，定义预估的初值。更好的方法是先对函数进行数学变换化成易解的形式再进行拟合。这个目前还没做。
总结：第一次用python数学建模，和matlab比起来其实差不多。因为有很多包可以用，有的包感觉基本就是把matlab的功能移植过来了，连用法和参数都差不多。但如果只是数学建模感觉还是matlab好用一点，至少省了调包的步骤。
原文：https://www.cnblogs.com/raybu/p/12554576.html
