【第9回Python流体の数値計算】Numpyでスライスを使って計算を速くする！

こんにちは(@t_kun_kamakiri)。

本日は、「Pythonでスライスによる要素の参照」というのをやります。

これにより、メインプログラムでfor文で計算を回すよりも、計算速度が速くなるので計算時間も比較しながら解析をしたいと思います。

本件の基本的な内容はこちらのサイトにそってやっていきます。

この記事ではこんな人を対象にしています。

こんな人が対象

Pythonを使い始めたけどどう使うかわからない
流体の数値計算をはじめて勉強する人

本記事シリーズの最終目標は、ナビエストークス方程式をPythonで実装することですが、いきなりナビエストークス方程式を実装するのは難しいので各項の意味を確認しながら進めていきたいと思います。

今回の内容はこちら

スライスの基本
for文とスライスによる計算時間の計測の比較

カマキリ

本記事の内容を一言でいうと「スライス」の方が計算速いよ！ってことです。

前回の記事までで、「１次元の移流方程式」「1次元の拡散方程式」「1次元バーガース方程式」をPythonで実装するというのをやりましたが、今回はPythonでの計算速度を早くする方法のひとつであるスライスについて慣れておきたいと思います。

Pythonの使用環境は、Google Colaboratoryを使っています。

スライスの基本
Numpy配列でのスライスの例
for文を使った場合とスライスを使った場合の計算時間を計測
スライスと使うときはndarray配列で使う
1. ndarray配列は要素ごとの演算
2. リスト型は要素の結合
２次元の移流方程式の数値計算
1. for文を使う
2. スライスを使う
まとめ
Pythonの完全初心者は書籍で学ぶとよい
流体の勉強をしたい方
数値計算を勉強したい方

スライスの基本

スライスは、リスト型だけではなく文字列、タプルにも使うことができます。

【ポイント】
スライスは「リスト」「文字列」「タプル」に使うことができる。

スライスに基本は例えば以下のように使います。

以下のようにリストを用意します。

a = [0,1,2,3,4,5,6,7,8,9,10]

a

a = [0,1,2,3,4,5,6,7,8,9,10]

【結果】

[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

1	[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

普通のリストですね。

これをスライスを使って要素を抽出します。

それにはいくつか書き方があるのですが、基本的な書き方は[start:stop:step]です。

例えば以下のように書くと、

a[1:8:3]

a[1:8:3]

【結果】

[1, 4, 7]

[1, 4, 7]

このようになります。
これは次のようになっていると理解してください。

スライスの基本的な書き方は上記となりますが、その他にも色々な書き方があるのでちょっと見ておきましょうかね(^^♪

すべての要素を取り出すときは以下のように書きます。

a[:]

a[:]

【結果】

[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

1	[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

「Index：1～３」の値を抽出したい場合はstepを指定しなければ、stepは１となるので以下のよう書けばOKです。

a[0:3]

a[0:3]

【結果】

[0, 1, 2]

[0, 1, 2]

また、次のように最後の２つ目まですべて抽出したい場合は、マイナスをつけることで抽出できます。

a[:-2]

a[:-2]

【結果】

[0, 1, 2, 3, 4, 5, 6, 7, 8]

1	[0, 1, 2, 3, 4, 5, 6, 7, 8]

基本的にはこのようにスライスを使って要素を抽出ることができます。

本記事ではリスト型を主に近いNumPy配列オブジェクトのndarrayでのスライスを使った解説がメインになります。

カマキリ

Pythonのコードを書きながら進めましょう！

Numpy配列でのスライスの例

まずは、以下のようなシンプルな例を考えてみましょう。

$\begin{align*} u^{n+1}_{i}=u^{n}_{i}-u^{n}_{i-1}\tag{1} \end{align*}$

$u^{n}=[0,1,2,3,4,5]$ なので、これを(1)式に従って計算する場合、今までは以下のようにfor文を使って計算していましたよね。

import numpy as np

u = np.array((0, 1, 2, 3, 4, 5))

for i in range(1, len(u)):
    print(u[i] - u[i-1])

import numpy as np

u = np.array((0, 1, 2, 3, 4, 5))

for i in range(1, len(u)):

print(u[i] - u[i-1])

【結果】

これは毎回、uという配列の要素をひとつひとつ取り出して、 u[i] – u[i-1]という計算を行っているので時間がかかるんですよね。

しかし、スライスというのを使えば、一回のコマンドで同様の計算ができて、かつ計算時間が短くて済みます。

カマキリ

たった1行です！

u[1:] - u[0:-1]

1	u[1:] - u[0:-1]

以上です。

これで先ほどのfor文と同じ結果になります。

これはndarray配列の場合は、各要素ごとに計算をすることになるので以下のように計算がされています。

for文を使った場合とスライスを使った場合の計算時間を計測

Pythonは上で見たようなスライスを使った方法の方が、for文を使うよりも計算が速くなります。

試しにfor文とスライスでの計算時間を計測して比較してみましょう。

コードの先頭に「%time」や「%%time」と書くことで、計算時間の計測してくれます。

%time：1回の計算時間を算出
%%time：何回か計算してその計算時間の平均を算出

では、まずはfor文から・・・

%%timeit
import numpy as np

u = np.array((0, 1, 2, 3, 4, 5))
for i in range(1, len(u)):
  u[i] - u[i-1]

%%timeit

import numpy as np

u = np.array((0, 1, 2, 3, 4, 5))

for i in range(1, len(u)):

u[i] - u[i-1]

【結果】

The slowest run took 7.71 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 3.63 µs per loop

1 2	The slowest run took 7.71 times longer than the fastest. This could mean that an intermediate result is being cached. 100000 loops, best of 3: 3.63 µs per loop

100000回の中でベストは「3.63μs」という結果ですね。

では、次にスライスを使ったコードを見てみましょう。

%%timeit
import numpy as np

u = np.array((0, 1, 2, 3, 4, 5))

u[1:] - u[0:-1]

%%timeit

import numpy as np

u = np.array((0, 1, 2, 3, 4, 5))

u[1:] - u[0:-1]

【結果】

The slowest run took 15.39 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 2.42 µs per loop

1 2	The slowest run took 15.39 times longer than the fastest. This could mean that an intermediate result is being cached. 100000 loops, best of 3: 2.42 µs per loop

100000回の中でベストは「2.42μs」という結果ですね。

やはり、スライスを使った方が計算時間が速いという結果になりましたね。

このように、Pythonはインタプリタ型の言語であるためコンパイラによる最適化がなく、コードの書き方の工夫で計算時間が速くなります。

Numpy自体が数値計算のためのライブラリであるため計算速度向上に貢献しているわけですが、for文を使ったごり押しの計算よりスマートにスライスと使えば、記述方法も計算時間も速くなるというメリットがあります。

少し余談になりますが、for文でndarray配列を作るときは、できるだけappendは使わない方が計算が速かったりもします。

例えば以下の2つのコードは同じ結果になりますが、計算時間が違います。

例１）appendを使って要素を格納

%%timeit
import numpy as np

du = []
u = np.array((0, 1, 2, 3, 4, 5))
for i in range(1, len(u)):
  du.append(u[i] - u[i-1])

%%timeit

import numpy as np

du = []

u = np.array((0, 1, 2, 3, 4, 5))

for i in range(1, len(u)):

du.append(u[i] - u[i-1])

【結果】

The slowest run took 48.42 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 4.24 µs per loop

1 2	The slowest run took 48.42 times longer than the fastest. This could mean that an intermediate result is being cached. 100000 loops, best of 3: 4.24 µs per loop

これは、空のリストファイル（du = []）を用意しておいて、appendを使って要素をリストに格納しています。
でも、この書き方はappendを毎回呼び出しているため計算時間のロスが生じているのですよね。

例2）appendを使わない

%%timeit
import numpy as np

du = []
u = np.array((0, 1, 2, 3, 4, 5))

du = np.array([i for i in u])

%%timeit

import numpy as np

du = []

u = np.array((0, 1, 2, 3, 4, 5))

du = np.array([i for i in u])

【結果】

The slowest run took 8.96 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 4.02 µs per loop

1 2	The slowest run took 8.96 times longer than the fastest. This could mean that an intermediate result is being cached. 100000 loops, best of 3: 4.02 µs per loop

わずかですが、appendを使うときよりも計算時間が速くなっていますよね。

これくらいの計算規模だと気にしない程度の計算速度ですが、大規模な計算になったときに計算速度の最適化は考えてコードを組んだ方が良いので、色々試行錯誤をして「こうすると速いよ」っていう書き方は覚えておく方がいいでしょう。

以上、余談でした！

スライスと使うときはndarray配列で使う

Pythonの勉強し始めのころは、要素ごとの演算のつもりがなぜか要素の結合になってしまったり、ndarray配列とリスト配列を全然区別できていませんでした。

カマキリ

当時の自分は雑魚だったってことです。

【ポイント】

ndarray配列は要素ごとの演算
リスト型は要素の結合

また、少し余談になりますが大事なところなので、解説しておこうと思います。

ndarray配列は要素ごとの演算

まずは、ndarray配列の生成方法を復習します。

ndarrayオブジェクトの生成方法は以下の3つを覚えておけば良いでしょう(^^)/
※はじめ、慣れないうちはtypeでデータの型を確認しながら進めた方が良いです。

例１）arrayを使う

import numpy as np

list = [0, 1, 2, 3, 4, 5]

u = np.array(list)

print(u)
print(type(u))

import numpy as np

list = [0, 1, 2, 3, 4, 5]

u = np.array(list)

print(u)

print(type(u))

【結果】

[0 1 2 3 4 5]
<class 'numpy.ndarray'>

1 2	[0 1 2 3 4 5] <class 'numpy.ndarray'>

例２）linspace(始め：終わり：要素数)を使う

import numpy as np

u = np.linspace(0,10,11)

print(u)
print(type(u))

import numpy as np

u = np.linspace(0,10,11)

print(u)

print(type(u))

【結果】

[ 0.  1.  2.  3.  4.  5.  6.  7.  8.  9. 10.]
<class 'numpy.ndarray'>

1 2	[ 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.] <class 'numpy.ndarray'>

例３）arange(始め：終わり：分割幅)を使う

import numpy as np

u = np.arange(0,10,1)

print(u)
print(type(u))

import numpy as np

u = np.arange(0,10,1)

print(u)

print(type(u))

【結果】

[0 1 2 3 4 5 6 7 8 9]
<class 'numpy.ndarray'>

1 2	[0 1 2 3 4 5 6 7 8 9] <class 'numpy.ndarray'>

さらに詳しく知りたい方は、以下の記事も参考にしてください。

Numpyについての記事はこちら

【第2回Python流体の数値計算】Numpyの使い方を確認する。

ndarrayオブジェクトを使えば、以下のように要素ごとの演算を行ってくれます。

import numpy as np

u1 = np.linspace(0,10,5)
u2 = np.linspace(0,20,5)

u3 = u1 + u2

print('u1 = ', u1)
print('u2 = ', u2)
print('u3 = ', u3)

import numpy as np

u1 = np.linspace(0,10,5)

u2 = np.linspace(0,20,5)

u3 = u1 + u2

print('u1 = ', u1)

print('u2 = ', u2)

print('u3 = ', u3)

【結果】

u1 =  [ 0.   2.5  5.   7.5 10. ]
u2 =  [ 0.  5. 10. 15. 20.]
u3 =  [ 0.   7.5 15.  22.5 30. ]

u1 = [ 0. 2.5 5. 7.5 10. ]

u2 = [ 0. 5. 10. 15. 20.]

u3 = [ 0. 7.5 15. 22.5 30. ]

めちゃくちゃ直感的に計算できて便利ですよね。

しかし、注意点として配列の要素数同士の演算を行うとエラーがおこるので注意してください。

例えば、以下のコードはエラーが生じます。

import numpy as np

u1 = np.linspace(0,10,5)
u2 = np.linspace(0,20,3)

u3 = u1 + u2

print('u1 = ', u1)
print('u2 = ', u2)
print('u3 = ', u3)

import numpy as np

u1 = np.linspace(0,10,5)

u2 = np.linspace(0,20,3)

u3 = u1 + u2

print('u1 = ', u1)

print('u2 = ', u2)

print('u3 = ', u3)

【結果】

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-31-d9d3a3b2b5ca> in <module>()
      4 u2 = np.linspace(0,20,3)
      5 
----> 6 u3 = u1 + u2
      7 
      8 print('u1 = ', u1)

ValueError: operands could not be broadcast together with shapes (5,) (3,)

---------------------------------------------------------------------------

ValueError Traceback (most recent call last)

<ipython-input-31-d9d3a3b2b5ca> in <module>()

4 u2 = np.linspace(0,20,3)

----> 6 u3 = u1 + u2

8 print('u1 = ', u1)

ValueError: operands could not be broadcast together with shapes (5,) (3,)

最後の行を見ると、要素の大きさが違うってエラーが書かれていますよね。

失敗した原因は「u1 = np.linspace(0,10,5)」は要素数が5個に対して、「u2 = np.linspace(0,20,3)」は要素数が3個だからです。
エラーの内容は要素数の異なるもの同士を足し算しようとしたからですね。

こんな感じで、ndarrayの演算は要素の数にお気を付けくださいってことです！

リスト型は要素の結合

次に、リストを用意します。

list = [0, 1, 2, 3, 4, 5]

print(list)
print(type(list))

list = [0, 1, 2, 3, 4, 5]

print(list)

print(type(list))

【結果】

[0, 1, 2, 3, 4, 5]
<class 'list'>

1 2	[0, 1, 2, 3, 4, 5] <class 'list'>

ここで注意なのですが、リスト型同士の足し算はリストの要素の結合になってしまいます。

list_1 = [0, 1, 2, 3, 4, 5]
list_2 = [0, 1, 2, 3, 4, 5]

list_3 = list_1 + list_2
print(list_3)
print(type(list_3))

list_1 = [0, 1, 2, 3, 4, 5]

list_2 = [0, 1, 2, 3, 4, 5]

list_3 = list_1 + list_2

print(list_3)

print(type(list_3))

【結果】

[0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5]
<class 'list'>

1 2	[0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5] <class 'list'>

なので、リスト同士の足し算は、各リストの要素の数は違ってもエラーになりません。

list_1 = [0, 1, 2]
list_2 = [0, 1, 2, 3, 4, 5]

list_3 = list_1 + list_2
print(list_3)
print(type(list_3))

list_1 = [0, 1, 2]

list_2 = [0, 1, 2, 3, 4, 5]

list_3 = list_1 + list_2

print(list_3)

print(type(list_3))

【結果】

[0, 1, 2, 0, 1, 2, 3, 4, 5]
<class 'list'>

1 2	[0, 1, 2, 0, 1, 2, 3, 4, 5] <class 'list'>

しかし、ndarrayのノリで以下の計算をリストのまま行うと、要素の結合になってしまうので注意が必要です。

リストのまま演算をしてしまった例

import numpy as np

u = [0, 1, 2, 3, 4, 5]
list = u[1:] + u[0:-1]
print(list)

import numpy as np

u = [0, 1, 2, 3, 4, 5]

list = u[1:] + u[0:-1]

print(list)

【結果】

[1, 2, 3, 4, 5, 0, 1, 2, 3, 4]

1	[1, 2, 3, 4, 5, 0, 1, 2, 3, 4]

こんな感じで「思ってた計算結果と違う！」ってなるのでここはよく理解しておきましょう。

ndarray配列で演算

import numpy as np

u =np.array( [0, 1, 2, 3, 4, 5])
list = u[1:] + u[0:-1]
print(list)

import numpy as np

u =np.array( [0, 1, 2, 3, 4, 5])

list = u[1:] + u[0:-1]

print(list)

【結果】

[1 3 5 7 9]

1	[1 3 5 7 9]

やりたいのはこんな感じで「要素ごとの演算」ですよね。

説明が長くなりましたが、初学者だった3年くらい前の自分はこのあたりを理解できていなかったので、詳しくまとめておきました。

では、今度は、スライスを使って実際に2次元の移流方程式を数値計算で解いてみます。

２次元の移流方程式の数値計算

１次元の移流方程式については、以下の記事でまとめてあります。

前回の記事はこちら

で、今回は２次元の移流方程式に拡張して数値計算を実装します。

$x$ 方向、 $y$ 方向の流れに対して流速を $u(x,y,t$ と $v(x,y,t$ と書くと、解くべき方程式は以下になります。

$\begin{align*} \frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=0\\ \frac{\partial v}{\partial t}+c\frac{\partial v}{\partial x}=0\tag{1} \end{align*}$

これを離散化して数値計算を行います。移流項の離散化は前回の記事でも同じで後退差分で行います。

$\begin{align*} u^{n+1}_{i,j}=u^{n}_{i,j}-\frac{c\Delta t}{\Delta x}\big(u^{n}_{i,j}-u^{n}_{i-1,j}\big)-\frac{c\Delta t}{\Delta y}\big(u^{n}_{i,j}-u^{n}_{i,j-1}\big)\\ v^{n+1}_{i,j}=v^{n}_{i,j}-\frac{c\Delta t}{\Delta y}\big(v^{n}_{i,j}-v^{n}_{i-1,j}\big)-\frac{c\Delta t}{\Delta y}\big(v^{n}_{i,j}-v^{n}_{i,j-1}\big)\tag{2} \end{align*}$

初期状態を設定して、for文とスライスを使って解く場合の計算時間を比較してみようと思います。

【手順】

初期状態の設定
（例１）for文で解く
（例２）スライスで解く

では、まず正気状態を設定します。

import numpy as np
nx = 81
ny = 81
nt = 100
c = 1
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
sigma = .2
dt = sigma * dx

x = np.linspace(0, 2, nx)
y = np.linspace(0, 2, ny)

u = np.ones((ny, nx)) ##create a 1xn vector of 1's
un = np.ones((ny, nx)) 

###Assign initial conditions

u[int(.5 / dy): int(1 / dy + 1), int(.5 / dx):int(1 / dx + 1)] = 2

import numpy as np

nx = 81

ny = 81

nt = 100

c = 1

dx = 2 / (nx - 1)

dy = 2 / (ny - 1)

sigma = .2

dt = sigma * dx

x = np.linspace(0, 2, nx)

y = np.linspace(0, 2, ny)

u = np.ones((ny, nx)) ##create a 1xn vector of 1's

un = np.ones((ny, nx))

###Assign initial conditions

u[int(.5 / dy): int(1 / dy + 1), int(.5 / dx):int(1 / dx + 1)] = 2

どんな初期状態になっているのか確認をしておきましょう。

import matplotlib.pyplot as plt

plt.plot(x,u[:,int(.5 / dy)])

import matplotlib.pyplot as plt

plt.plot(x,u[:,int(.5 / dy)])

これで初期状態を作ることができました。

for文を使う

次に普通にfor文を使って時間に対して0～nt+1ステップの計算をさせて計算時間を測定してみます。
0～nt+1ステップの1ステップごとにx−y空間のfor文を2回使うことになりますね。
考えただけでも計算時間がかかりそうなアルゴリズムです・・・

%%timeit
u = np.ones((ny, nx))
u[int(.5 / dy): int(1 / dy + 1), int(.5 / dx):int(1 / dx + 1)] = 2

for n in range(nt + 1): ##loop across number of time steps
    un = u.copy()
    row, col = u.shape
    for j in range(1, row):
        for i in range(1, col):
            u[j, i] = (un[j, i] - (c * dt / dx * 
                                  (un[j, i] - un[j, i - 1])) - 
                                  (c * dt / dy * 
                                   (un[j, i] - un[j - 1, i])))
            u[0, :] = 1
            u[-1, :] = 1
            u[:, 0] = 1
            u[:, -1] = 1

%%timeit

u = np.ones((ny, nx))

u[int(.5 / dy): int(1 / dy + 1), int(.5 / dx):int(1 / dx + 1)] = 2

for n in range(nt + 1): ##loop across number of time steps

un = u.copy()

row, col = u.shape

for j in range(1, row):

for i in range(1, col):

u[j, i] = (un[j, i] - (c * dt / dx *

(un[j, i] - un[j, i - 1])) -

(c * dt / dy *

(un[j, i] - un[j - 1, i])))

u[0, :] = 1

u[-1, :] = 1

u[:, 0] = 1

u[:, -1] = 1

※計算時間の計測のためなのて、vについては今回は解いていません（２次元の数値計算については詳しくは次回の記事で解説をします）

1 loop, best of 3: 2.96 s per loop

1	1 loop, best of 3: 2.96 s per loop

計算時間「2.96 s」でした。約３秒かかるということですね。

スライスを使う

次に、スライスを使って計算を実行してみます。

%%timeit
u = np.ones((ny, nx))
u[int(.5 / dy): int(1 / dy + 1), int(.5 / dx):int(1 / dx + 1)] = 2

for n in range(nt + 1): ##loop across number of time steps
    un = u.copy()
    u[1:, 1:] = (un[1:, 1:] - (c * dt / dx * (un[1:, 1:] - un[1:, 0:-1])) -
                              (c * dt / dy * (un[1:, 1:] - un[0:-1, 1:])))
    u[0, :] = 1
    u[-1, :] = 1
    u[:, 0] = 1
    u[:, -1] = 1