這學期選了數學建模課，因為上學期學了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好用一點，至少省了調包的步驟。