这学期选了数学建模课，因为上学期学了MATLAB，这学期尝试使用Python完成数学建模。
Python的基本语法其实很简单
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)View Code
这里借鉴了一下网上的代码。主函数是全部自己写的。遇到的问题是：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好用一点，至少省了调包的步骤。
