2. Numpy array (配列) のコツ.

python では listという概念がありますが,c言語やfotranで習う配列とはかなり異なる概念です. listの使い方については python user会 を参照して下さい.

2.1. 基本的な事

数値計算をする上でlistを配列の様に使う事ができますが,

>>> l = [1,2,3,4,5]
>>> g = [6,7,8,9,10]
>>> for i in range(5):
...     a = l[i] + g[i]
...     print(a)
...
7
9
11
13
15

listを配列の様に使うと非常に計算が遅くなります.そこで一般にはnumpyの配列 arrayがよく使われます. numpy.array 正確には numpy.ndarrayは (固定長の)N次元の配列を使う為のインターフェイスで,Cの配列以上に高機能です. 正確な事が知りたければ numpy reference を参照する事をお勧めします.

基本的な使い方は日本語でも検索すればいくらでも出てくると思いますが,幾つか簡単に示します. 上の例をnumpy.arrayをつかって書くと

>>> import numpy as np
>>> a = np.array([1,2,3,4,5])
>>> b = np.array([6,7,8,9,10])
>>> for i in range(5):
...     print(a[i] + b[i])
...
7
9
11
13
15

ですが,各要素の足し算であれば

>>> a + b
array([ 7,  9, 11, 13, 15])

で十分でしょう.むしろfor文を使って要素を一つづつ計算すると計算時間が長くなります.

>>> a = np.random.random(10**7)
>>> b = np.random.random(10**7)
>>> %time for i in range(10**6): a[i] + b[i]
CPU times: user 601 ms, sys: 741 µs, total: 602 ms
Wall time: 601 ms
>>> %time a + b
CPU times: user 21.6 ms, sys: 21.9 ms, total: 43.5 ms
Wall time: 47.7 ms
array([ 0.97029472,  1.25966728,  0.75212228, ...,  1.63341271,
        0.51066578,  1.42720649])

ここで%timeは ipythonの magic function です.

2.2. 数値計算でよく使うテクニック

私がよく使うテクニックを記載しておきます.一般的かどうかは知りません.

>>> np.array([0]*5)
array([0, 0, 0, 0, 0])
>>> np.zeros(5,dtype=np.int)
array([0, 0, 0, 0, 0])

両者は同じものだが,前者は任意の要素をarrayに格納できる.

arrayはpythonの list内包表記 が可能です.

>>> np.array([i for i in range(10)])
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> np.array([i if i%2 ==0 else 0 for i in range(10)])
array([0, 0, 2, 0, 4, 0, 6, 0, 8, 0])
>>> np.array([[2]*5 for i in range(5)])
array([[2, 2, 2, 2, 2],
       [2, 2, 2, 2, 2],
       [2, 2, 2, 2, 2],
       [2, 2, 2, 2, 2],
       [2, 2, 2, 2, 2]])

常微分方程式や偏微分方程式の計算の際に必要となる差分計算.例えば

\[\begin{split}\frac{\partial u}{\partial x}&\to \frac{u(x+\Delta x) - u(x)}{\Delta x}\\ \frac{\partial^2 u}{\partial x^2}&\to \frac{u(x+\Delta x) - 2u(x) + u(x-\Delta x)}{(\Delta x)^2}\end{split}\]

のような差分が良く現れる.素朴に計算すると(以下の例では分母を無視する)

>>> u = np.array([np.random.randint(10) for i in range(10)])
>>> print(u)
[6 1 5 5 1 3 6 8 5 8]
>>> res1 = np.array([0.0]*9)
>>> res2 = np.array([0.0]*8)
>>> for i in range(0,len(u)-1):
...     res1[i] = u[i+1] - u[i]
...
>>> print(res1)
[-5.  4.  0. -4.  2.  3.  2. -3.  3.]
>>> for i in range(1,len(u)-1):
...     res2[i-1] = u[i+1] - 2*u[i] + u[i-1]
...
>>> print(res2)
[ 9. -4. -4.  6.  1. -1. -5.  6.]

であるが,これはfor文で記述されているので配列の次元が増えると計算に時間が掛かる. numpyには差分( numpy.diff )関数が容易されており

>>> np.diff(u)
array([-5,  4,  0, -4,  2,  3,  2, -3,  3])
>>> np.diff(u,n=2)
array([ 9, -4, -4,  6,  1, -1, -5,  6])

とする方が高速であり簡潔である.1次元配列の場合はdiff関数でも良いが,次の方法が多次元配列を扱う上で 拡張性が容易である.

>>> u[1:] - u[0:-1]
array([-5,  4,  0, -4,  2,  3,  2, -3,  3])
>>> u[2:] -2*u[1:-1] + u[:-2]
array([ 9, -4, -4,  6,  1, -1, -5,  6])

実際 \(u(x,y)\) の場合,\(\nabla u, \nabla^2 u\) の陽的差分

\[\begin{split}\frac{\partial u}{\partial x} + \frac{\partial u}{\partial y}&\to \frac{u(x+\Delta x,y) - u(x,y)}{\Delta x}+ \frac{u(x,y+\Delta y) - u(x,y)}{\Delta y}\\\\ \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} &\to \frac{u(x+\Delta x, y) - 2u(x,y) + u(x-\Delta x,y)}{(\Delta x)^2} + \frac{u(x,y+\Delta y) - 2u(x,y) + u(x,y-\Delta y)}{(\Delta y)^2}\end{split}\]

を素朴に計算すると

u = np.random.randint(10,size=(5,5))
>>> res1 = np.array([[0.0]*4 for i in range(4)])
>>> res2 = np.array([[0.0]*3 for i in range(3)])
>>> u
array([[0, 5, 1, 5, 9],
       [7, 2, 7, 7, 6],
       [6, 6, 3, 1, 1],
       [3, 5, 9, 8, 0],
       [6, 5, 5, 4, 0]])
>>> for i in range(0, u.shape[0]-1):
        for j in range(0, u.shape[1]-1):
            res1[i,j] = u[i+1,j] - u[i,j] + u[i,j+1] - u[i,j]
...
>>> res1
array([[ 12.,  -7.,  10.,   6.],
       [ -6.,   9.,  -4.,  -7.],
       [ -3.,  -4.,   4.,   7.],
       [  5.,   4.,  -5., -12.]])
>>> for i in range(1,u.shape[0]-1):
        for j in range(1,u.shape[1]-1):
            res2[i-1,j-1] = u[i+1,j] - 2*u[i,j] + u[i-1,j] + u[i,j+1] - 2*u[i,j]+ u[i,j-1]
...
>>> res2
array([[ 17., -15.,  -9.],
       [ -8.,  11.,  15.],
       [  3., -15., -18.]])

であるが,配列の次元だけfor分が増えるのでpythonでは絶望的に遅くなる.故に次の方法が最も効率的だと思われる.

>>> u[1:,0:-1] - u[0:-1,0:-1] + u[0:-1,1:] -u[0:-1,0:-1]
array([[ 12,  -7,  10,   6],
       [ -6,   9,  -4,  -7],
       [ -3,  -4,   4,   7],
       [  5,   4,  -5, -12]])
>>> (u[2:,1:-1] - 2*u[1:-1,1:-1] + u[0:-2,1:-1] +
 u[1:-1,2:] - 2*u[1:-1,1:-1] + u[1:-1,0:-2])
>>>
array([[ 17, -15,  -9],
       [ -8,  11,  15],
       [  3, -15, -18]])

これなら多次元への拡張が容易であり,速度も問題ない.diffを使って同等のこともできるかもしれないが,私は方法を知らない.