數學建模三劍客MSN

前言

不論是不是巴薩的球迷，只要你喜歡足球，就必定據說過梅西(Messi)、蘇亞雷斯(Suarez)和內馬爾(Neymar)這個MSN組合。在衆多的數學建模輔助工具中，也有一個犀利無比的MSN組合，他們就是python麾下大名鼎鼎的 Matplotlib + Scipy + Numpy三劍客。

本文是我整理的MSN學習筆記，有些理解可能比較膚淺，甚至是錯誤的。若是所以誤導了某位看官，在工做中形成重大失誤或損失，我頂多只能賠償一頓飯——還得是咱們樓下的十元盒飯。特此聲明。

文中代碼均從個人這臺時不時出點問題、鬧個情緒的Yoga 3 pro上覆制而來，這意味着全部的代碼都可在下面的運行環境中順利運行：

• pyhton 2.7.8
• numpy 1.11.1
• scipy 0.16.1
• matplotlib 1.5.1

三劍客之Numpy

numpy是一個開源的python科學計算庫，包含了不少實用的數學函數，涵蓋線性代數、傅里葉變換和隨機數生成等功能。最初的numpy實際上是scipy的一部分，後來才從scipy中分離出來。

numpy不是python的標準庫，須要單獨安裝。假定你的運行環境已經安裝了python包管理工具pip，numpy的安裝就很是簡單：

pip install numpy

數組對象

ndarray是多維數組對象，也是numpy最核心的對象。在numpy中，數組的維度(dimensions)叫作軸(axes)，軸的個數叫作秩(rank)。一般，一個numpy數組的全部元素都是同一種類型的數據，而這些數據的存儲和數組的形式無關。

下面的例子，建立了一個三維的數組（在導入numpy時，通常都簡寫成np）。

import numpy asnp

a= np.array([[1,2,3],[4,5,6],[7,8,9]])

數據類型

numpy支持的數據類型主要有布爾型(bool)、整型(integrate)、浮點型(float)和複數型(complex)，每一種數據類型根據佔用內存的字節數又分爲多個不一樣的子類型。常見的數據類型見下表。

建立數組

一般，咱們用np.array()建立數組。若是僅僅是建立一維數組，也可使用np.arange()或者np.linspace()的方法。np.zeros()、np.ones()、np.eye()則能夠構造特殊的數據。np.random.randint()和np.random.random()則能夠構造隨機數數組。

構造複雜數組

不少時候，咱們須要從簡單的數據結構，構造出複雜的數組。例如，用一維的數據生成二維格點。

重複數組: tile

一維數組網格化: meshgrid

指定範圍和分割方式的網格化: mgrid

上面的例子中用到了虛數。構造虛數的方法以下：

>>> complex(2,5)

(2+5j)

數組的屬性

numpy的數組對象除了一些常規的屬性外，也有幾個相似轉置、扁平迭代器等看起來更像是方法的屬性。扁平迭代器也許是遍歷多維數組的一個簡明方法，下面的代碼給出了一個例子。

改變數組維度

numpy數組的存儲順序和數組的維度是不相干的，所以改變數組的維度是很是便捷的操做，除resize()外，這一類操做不會改變所操做的數組自己的存儲順序。

索引和切片

對於一維數組的索引和切片，numpy和python的list同樣，甚至更靈活。

假設有一棟2層樓，每層樓內的房間都是3排4列，那咱們能夠用一個三維數組來保存每一個房間的居住人數（固然，也能夠是房間面積等其餘數值信息）。

數組合並

數組合併除了下面介紹的水平合併、垂直合併、深度合併外，還有行合併、列合併，以及concatenate()等方式。假如你比我還懶，那就只瞭解前三種方法吧，足夠用了。

數組拆分

拆分是合併的逆過程，概念是同樣的，但稍微有一點不一樣：

數組運算

數組和常數的四則運算，是數組的每個元素分別和常數運算；數組和數組的四則運算則是兩個數組對應元素的運算（兩個數組有相同的shape，不然拋出異常）。

特別提示：若是想對數組內符合特定條件的元素作特殊處理，下面的代碼也許有用。

數組方法和經常使用函數

數組對象自己提供了計算算數平均值、求最大最小值等內置方法，numpy也提供了不少實用的函數。爲了縮減篇幅，下面的代碼僅以一維數組爲例，展現了這些方法和函數用法。事實上，大多數狀況下這些方法和函數對於多維數組一樣有效，只有少數例外，好比compress函數。

矩陣對象

matrix是矩陣對象，繼承自ndarray類型，所以含有ndarray的全部數據屬性和方法。不過，當你把矩陣對象當數組操做時，須要注意如下幾點：

• matrix對象老是二維的，即便是展平（ravel函數）操做或是成員選擇，返回值也是二維的
• matrix對象和ndarray對象混合的運算老是返回matrix對象

建立矩陣

matrix對象可使用一個Matlab風格的字符串來建立（以空格分隔列，以分號分隔行的字符串），也能夠用數組來建立。

矩陣的特有屬性

矩陣有幾個特有的屬性使得計算更加容易，這些屬性有：

矩陣乘法

對ndarray對象而言，星號是按元素相乘，dot()函數則看成矩陣相乘。對於matrix對象來講，星號和dot()函數都是矩陣相乘。特別的，對於一維數組，dot()函數實現的是向量點乘（結果是標量），但星號實現的卻不是差乘。

線性代數模塊

numpy.linalg 是numpy的線性代數模塊，能夠用來解決逆矩陣、特徵值、線性方程組以及行列式等問題。

計算逆矩陣

儘管matrix對象自己有逆矩陣的屬性，但用numpy.linalg模塊求解矩陣的逆，也是很是簡單的。

計算行列式

如何計算行列式，我早已經不記得了，但手工計算行列式的痛苦，我依然記憶猶新。如今好了，你在手機上均可以用numpy輕鬆搞定（前提是你的手機上安裝了python + numpy）。

m= np.mat('0 1 2; 1 0 3; 4 -3 8')

np.linalg.det(m)# 什麼？這就成了？

2.0

計算特徵值和特徵向量

截至目前，個人工做和特徵值、特徵向量還有沒任何關聯。記錄這一節，純粹是爲了我女兒，她正在讀數學專業。

求解線性方程組

有線性方程組以下：

x – 2y + z = 0

2y -8z = 8

-4x + 5y + 9z = -9

求解過程以下：

三劍客之Matplotlib

matplotlib 是python最著名的繪圖庫，它提供了一整套和Matlab類似的命令API，十分適合交互式地進行製圖。並且也能夠方便地將它做爲繪圖控件，嵌入GUI應用程序中。matplotlib 能夠繪製多種形式的圖形包括普通的線圖，直方圖，餅圖，散點圖以及偏差線圖等；能夠比較方便的定製圖形的各類屬性好比圖線的類型，顏色，粗細，字體的大小等；它可以很好地支持一部分 TeX 排版命令，能夠比較美觀地顯示圖形中的數學公式。

pyplot介紹

Matplotlib 包含了幾十個不一樣的模塊， 如 matlab、mathtext、finance、dates 等，而 pyplot 則是咱們最經常使用的繪圖模塊，這也是本文介紹的重點。

中文顯示問題的解決方案

有不少方法能夠解決此問題，但下面的方法恐怕是最簡單的解決方案了（我只在windows平臺上測試過，其餘平臺請看官自測）。

繪製最簡單的圖形

>>> import numpy asnp

>>> import matplotlib.pyplot asplt

>>> x= np.arange(0,2*np.pi,0.01)

>>> y= np.sin(x)

>>> plt.plot(x,y)

>>> plt.show()

設置標題、座標軸名稱、座標軸範圍

若是你在python的shell中運行下面的代碼，而shell的默認編碼又不是utf-8的話，中文可能仍然會顯示爲亂碼。你能夠嘗試着把 u’正弦曲線’ 寫成 ‘正弦曲線’.decode(‘gbk’)或者‘正弦曲線’.decode(‘utf-8’)

設置點和線的樣式、寬度、顏色

plt.plot函數的調用形式以下：

plot(x,y,color='green',linestyle='dashed',linewidth=1,marker='o',markerfacecolor='blue',markersize=6)

plot(x,y,c='g',ls='--',lw=1,marker='o',mfc='blue',ms=6)

1. color指定線的顏色，可簡寫爲「c」。顏色的選項爲：
• 藍色： ‘b’ (blue)
• 綠色： ‘g’ (green)
• 紅色： ‘r’ (red)
• 墨綠： ‘c’ (cyan)
• 洋紅： ‘m’ (magenta)
• 黃色： ‘y’ (yellow)
• 黑色： ‘k’ (black)
• 白色： ‘w’ (white)
• 灰度表示： e.g. 0.75 ([0,1]內任意浮點數)
• RGB表示法： e.g. ‘#2F4F4F’ 或 (0.18, 0.31, 0.31)
2. linestyle指定線型，可簡寫爲「ls」。線型的選項爲：
• 實線： ‘-’ (solid line)
• 虛線： ‘–’ (dashed line)
• 虛點線： ‘-.’ (dash-dot line)
• 點線： ‘:’ (dotted line)
• 無： 」或’ ‘或’None’
3. linewidth指定線寬，可簡寫爲「lw」。
4. marker描述數據點的形狀
• 點線： ‘.’
• 點線： ‘o’
• 加號： ‘+
• 叉號： ‘x’
• 上三角： ‘^’
• 上三角： ‘v’
5. markerfacecolor指定數據點標記的表面顏色，可 簡寫爲「 mfc」。
6. markersize指定數據點標記的大小，可 簡寫爲「 ms」。

文本標註和圖例

咱們分別使用不一樣的線型、顏色來繪製以十、e、2爲基的一組冪函數曲線，演示文本標註和圖例的使用。

在繪製圖例時，loc用於指定圖例的位置，可用的選項有：

• best
• upper right
• upper left
• lower left
• lower right

繪製多軸圖

在介紹如何將多幅子圖繪製在同一畫板的同時，順便演示如何繪製直線和矩形。咱們可使用subplot函數快速繪製有多個軸的圖表。subplot函數的調用形式以下：

subplot(numRows, numCols, plotNum)

subplot將整個繪圖區域等分爲numRows行 * numCols列個子區域，而後按照從左到右，從上到下的順序對每一個子區域進行編號，左上的子區域的編號爲1。若是numRows，numCols和plotNum這三個數都小於10的話，能夠把它們縮寫爲一個整數，例如subplot(323)和subplot(3,2,3)是相同的。subplot在plotNum指定的區域中建立一個軸對象。若是新建立的軸和以前建立的軸重疊的話，以前的軸將被刪除。

三劍客之Scipy

前面已經說過，最初的numpy實際上是scipy的一部分，後來才從scipy中分離出來。scipy函數庫在numpy庫的基礎上增長了衆多的數學、科學以及工程計算中經常使用的庫函數。例如線性代數、常微分方程數值求解、信號處理、圖像處理、稀疏矩陣等等。因爲其涉及的領域衆多，我之於scipy，就像盲人摸大象，只能是摸到哪兒算哪兒。

插值

一維插值和二維插值，是我最經常使用的scipy的功能之一，也是最容易上手的。

一維插值和樣條插值

下面的例子清楚地展現了線性插值和樣條插值以後的數據形態。

將原始數據以及線性插值和樣條插值以後的數據繪製在一塊兒，效果會比較明顯：

代碼以下：

特別說明：樣條插值附帶了不少默認參數，下面是簡單的說明。詳情請自行搜索。

scipy.interpolate.splrep(x,y,w=None,xb=None,xe=None,k=3,task=0,s=None,t=None,full_output=0,per=0,quiet=1)

# 參數s用來肯定平滑點數，一般是m-SQRT(2m),m是曲線點數。若是在插值中不須要平滑應該設定s=0。splrep()輸出的是一個3元素的元胞數組（t,c,k）,其中t是曲線點，c是計算出來的係數，k是樣條階數，一般是3階，但能夠經過k改變。

scipy.interpolate.splev(x,tck,der=0)

# 其中的der是進行樣條計算是須要實際計算到的階數，必須知足條件der<=k

擬合

在工做中，咱們經常須要在圖中描繪某些實際數據觀察的同時，使用一個曲線來擬合這些實際數據。所謂擬合，就是找出符合數據變化趨勢的曲線方程，或者直接繪製出擬合曲線。

使用numpy.polyfit擬合

下面這段代碼，基於Numpy模塊，能夠直接繪製出擬合曲線，但我沒法獲得曲線方程（儘管輸出了一堆曲線參數）。這是一個值得繼續深刻研究的問題。

3個擬合結果顯示在下圖中。

使用scipy.optimize.optimize.curve_fit擬合

scipy提供的擬合，貌似須要先肯定帶參數的曲線方程，而後由scipy求解方程，返回曲線參數。咱們仍是以上面的一組數據爲例使用scipy擬合曲線。

能夠看出，曲線近似正弦函數。構建函數y=a*sin(x*pi/6+b)+c，使用scipy的optimize.curve_fit函數求出a、b、c的值：

求解非線性方程（組）

在數學建模中，須要對一些稀奇古怪的方程（組）求解，Matlab天然是首選，但Matlab不是免費的，scipy則爲咱們提供了免費的午飯！scipy.optimize庫中的fsolve函數能夠用來對非線性方程（組）進行求解。它的基本調用形式以下：

fsolve(func, x0)

func(x)是計算方程組偏差的函數，它的參數x是一個矢量，表示方程組的各個未知數的一組可能解，func返回將x代入方程組以後獲得的偏差；x0爲未知數矢量的初始值。

咱們先來求解一個簡單的方程：sin(x)−cos(x)=0.2

>>> from scipy.optimize import fsolve

>>> import numpy asnp

>>> deff(A):

x= float(A[0])

return[np.sin(x)- np.cos(x)- 0.2]

>>> result= fsolve(f,[1])

array([0.92729522])

>>> print result

[0.92729522]

>>> printf(result)

[2.7977428707082197e-09]

哈哈，易如反掌！再來一個方程組：

4x2−2sin(yz)=0

5y+3=0

yz−1.5=0

圖像處理

在scipy.misc模塊中，有一個函數能夠載入Lena圖像——這副圖像是被用做圖像處理的經典示範圖像。我只是簡單展現一下在該圖像上的幾個操做。

1. 載入Lena圖像，並顯示灰度圖像
2. 應用中值濾波掃描信號的每個數據點，並替換爲相鄰數據點的中值
3. 旋轉圖像
4. 應用Prewitt濾波器（基於圖像強度的梯度計算）