Python數值計算:一 使用Pylab繪圖(1)

Pylab的使用

學習使用Python進行科學計算，然而很難找到簡單實用，又循序漸進的例子。正好手邊有一本《Matlab可視化大學物理學》，里面的例子非常清晰地解釋了Matlab在物理學中的應用。重新使用Python實現這些例子，學習了Python，又學習了科學計算思想，一舉兩得。

1.1Pylab簡介

pylab是Matplotlib和Ipython提供的一個模塊，提供了類似Matlab的語法，在啟動Ipython時可以使用--pylab啟動，它相當於導入了如下庫：

import numpy
import matplotlib
from matplotlib import pylab, mlab, pyplot
np = numpy
plt = pyplot
from IPython.display import display
from IPython.core.pylabtools import figsize, getfigs
from pylab import *
from numpy import *

1.2簡單應用

1.2.1畫圓

\(cos^2x+sin^2x = 1\)

n = 10
 #生成50個元素的行向量
phi = linspace(0, 2*pi)
 # 此時繪制x,y可以得到一個圓
x = cos(phi)
y = sin(phi)
for i in range(1,n):
    for j in range(1,i):
        #(繪圖 1)
        plot(x + 2 * j, y - 2 * j)

如果將繪圖1指令改為:

plot(x + 2*j -i, y - sqrt(3) * i)

可以繪制出交錯堆疊。

1.2.2不定循環（indefinite loop）與超越方程（transcentental equation）

超越方程在數學與物理中十分常見，然而並沒有確定的解析解法。需要使用數值解或圖形解。

在黑體輻射中，峰值波長與溫度成反比，比例系數為維恩常數（Wien's displacement law）。根據普朗克公式求維恩常數時，產生如下的超越方程：

\[x+5(e^{-x}-1)=0 \]

超越方程的解也應該是下面一條直線和指數函數曲線的交點：

\[y=x , y = 5[1-e^{(-x)}] \]

對於數值計算而言，可以使用牛頓法求解以上方程。

x0 = 1 #iteration method
xx = []
while True:
    x = 5 * (1 - exp( - x0))
    xx.append(x)
    if len(xx) > 1000:
        break
    if abs(x0 - x) < 0.0001:
        break
    x0 = x
subplot(211)#繪圖部分，首先繪制iteration method求得的解
plot(xx,"x-")
grid()
text(x0, x0, str(x0))
subplot(212)#然后繪制圖形解法的解
nx = arange(0,8,0.01)
plot(nx, nx, nx, 5*(1 - exp( -nx)), "--")
plot(x0, x0, "o")
text(x0, x0, str(x0))
grid()
show()