２次元の波¶
ライブラリの復習¶
Python は、Fortran、C、または C++ とは異なり、プログラミングを容易にすることを目的とした上位言語です。Pythonをインストールすると、豊富で用途の広い 標準ライブラリ が含まれており、個別のパッケージをダウンロードすることなく利用することできます。
ただし、標準ライブラリを使用することがいつも最適であるとは限りません。
必要に応じてライブラリをインストールして、効率よく作業を進めてください。
高度なライブラリを呼び出す方が良い場合もあります。以下は広く使用されているライブラリのリストです。
オペレーションシステム（OS）
os,systems時間
time,calendar数学
random,mathWebサービス
urllib,ftplib,xmlrpclib,jsonプロット
matplotlibグラフィカルユーザーインターフェイス
Tkinter,ttk科学計算
numpy,scipy機械学習
scikit-learnデータベース
dbm,gdmb,sqlite3
科学計算を行うNumpy¶
[7]:
import numpy as np # numpyを短い名前npにしてインポート
array = 5*np.ones(5) # サブ関数の呼び出し
print(array)
[5. 5. 5. 5. 5.]
Web ページのコンテンツを取得するurllib¶
[8]:
import urllib.request
req = urllib.request.Request('https://iikubo-lab.com')
response = urllib.request.urlopen(req)
the_page = response.read().decode()
#print(dir(the_page))
#print(the_page)
print(the_page.count('飯久保'))
16
使用するにはインポート¶
これらを使用するために必要なことは、インポートすることだけです。
import numpy as np # numpyを短い名前npにしてインポート
array = np.zeros(5) # サブ関数の呼び出し
または、
from numpy import zeros # zero関数を明示してライブラリから読み込み
array = zeros(5) # サブ関数の呼び出し
ない場合はインストール¶
ライブラリが存在しないというエラーメッセージが出力される場合には、anacondaプロンプトなどからインストールする必要があります。
$ conda install numpy
ライブラリは数が多くすべてを探索することは不可能です。目的にあったライブラリを見つけるには、google で検索してみてください。ここからは、いつもの numpy と matplotlib を使用して進めます。
クイズ¶
インポートしたライブラリで利用可能な関数を知るにはどうすればよいでしょうか。
Numpy, matplotlibの復習¶
この 2 つのライブラリはMathWorks社が開発している数値解析ソフトウェア MATLAB に非常に似ています。
NumPy は Python の科学計算の基本パッケージです。http://www.numpy.org/
以下のような特徴を持っています。 - N 次元配列オブジェクト - ブロードキャスト機能 - C/C++ と Fortran コードを統合するためのツール - 線形代数、フーリエ変換、乱数など
Matplotlib は高品質の図のための Python 2D プロット ライブラリです。https://matplotlib.org/
以下のような特徴を持っています。 - Python および IPython シェル - Jupyter notebook - Web アプリケーション サーバー - グラフィカルユーザーインターフェイス(GUI)ツールキット
簡単な使用例¶
配列の定義¶
[9]:
import numpy as np
array = np.linspace(0,1,10) #[Start, End, number of splits]#
print(array)
print(type(array))
[0. 0.11111111 0.22222222 0.33333333 0.44444444 0.55555556
0.66666667 0.77777778 0.88888889 1. ]
<class 'numpy.ndarray'>
他にもいろいろな方法があるので、以下の関数も試してみましょう。
[10]:
print(np.arange(0, 5, 0.1))
#print(np.array([1,2,3,4,5,6]))
#print(np.zeros([3,3]))
#print(np.ones([3,4]))
#print(np.eye(5))
[0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. 1.1 1.2 1.3 1.4 1.5 1.6 1.7
1.8 1.9 2. 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3. 3.1 3.2 3.3 3.4 3.5
3.6 3.7 3.8 3.9 4. 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9]
配列の操作¶
[11]:
import numpy as np
a = np.linspace(0,1,10) #[Start, End, number of splits]#
b = np.linspace(0,2,10) #[Start, End, number of splits]#
print(a)
print(b)
print(a+b)
print(1+b)
[0. 0.11111111 0.22222222 0.33333333 0.44444444 0.55555556
0.66666667 0.77777778 0.88888889 1. ]
[0. 0.22222222 0.44444444 0.66666667 0.88888889 1.11111111
1.33333333 1.55555556 1.77777778 2. ]
[0. 0.33333333 0.66666667 1. 1.33333333 1.66666667
2. 2.33333333 2.66666667 3. ]
[1. 1.22222222 1.44444444 1.66666667 1.88888889 2.11111111
2.33333333 2.55555556 2.77777778 3. ]
他にもいろいろな方法があるので、以下の関数も試してみましょう。
[12]:
a-b
#np.max(a)
#a+1.0
#a*b
#a_trans = a.transpose()
#np.dot(a_trans, b)
#dir(np.linalg)
#print(a_trans)
[12]:
array([ 0. , -0.11111111, -0.22222222, -0.33333333, -0.44444444,
-0.55555556, -0.66666667, -0.77777778, -0.88888889, -1. ])
サイン関数のグラフ¶
[13]:
import numpy as np
import matplotlib.pyplot as plt
x = np.arange(0, 10, 0.1)
y = np.sin(x)
plt.figure(figsize=(5,3))
plt.plot(x, y)
plt.show()
サイン関数のグラフ２¶
[14]:
pi = 3.1415926
plt.figure(figsize=(5,3))
plt.plot(x, y, '-.bd', label='sin(x)')
#plt.semilogx(x, y, 'dr')
#plt.loglog(x,y,'dr')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.xlim([0, 2*pi])
plt.ylim([-2, 2])
plt.show()
使い方がわからない場合にはヘルプを参照しましょう。
[15]:
help(plt.loglog)
Help on function loglog in module matplotlib.pyplot:
loglog(*args, **kwargs)
Make a plot with log scaling on both the x- and y-axis.
Call signatures::
loglog([x], y, [fmt], data=None, **kwargs)
loglog([x], y, [fmt], [x2], y2, [fmt2], ..., **kwargs)
This is just a thin wrapper around `.plot` which additionally changes
both the x-axis and the y-axis to log scaling. All the concepts and
parameters of plot can be used here as well.
The additional parameters *base*, *subs* and *nonpositive* control the
x/y-axis properties. They are just forwarded to `.Axes.set_xscale` and
`.Axes.set_yscale`. To use different properties on the x-axis and the
y-axis, use e.g.
``ax.set_xscale("log", base=10); ax.set_yscale("log", base=2)``.
Parameters
----------
base : float, default: 10
Base of the logarithm.
subs : sequence, optional
The location of the minor ticks. If *None*, reasonable locations
are automatically chosen depending on the number of decades in the
plot. See `.Axes.set_xscale`/`.Axes.set_yscale` for details.
nonpositive : {'mask', 'clip'}, default: 'clip'
Non-positive values can be masked as invalid, or clipped to a very
small positive number.
**kwargs
All parameters supported by `.plot`.
Returns
-------
list of `.Line2D`
Objects representing the plotted data.
その他の多くの使用例は、以下のサイトで確認することができます。
データからグラフを作る¶
初めからx値とy値が用意されておらず、一点一点計算してデータを追加するようなことがよくあります。そのような場合には、x座標とy座標の2つの空のリストを作成して、それらに値を追加します。ここではそのような形で、 \(sin(x)\) のグラフを作成する方法を示します。
[16]:
import matplotlib.pyplot as plt
from math import sin, cos
from numpy import linspace
#空のリストを作成
x = []
y1 = []
y2 = []
#リストに値を代入
for i in linspace(0,10,100):
x.append(i)
y1.append(sin(i))
y2.append(cos(i))
[17]:
#プロット
plt.plot(x, y1, 'r', label='sin(x)')
plt.plot(x, y2, 'o', label='cos(x)')
plt.legend()
plt.xlabel('x')
plt.ylabel('f(x)')
plt.show()
２次元の波の伝わり方¶
理工学分野では、2次元格子上のデータを扱う場合がしばしばあります。
固体物理：電荷や温度の変化、または固体表面上の原子堆積測定など
流体力学：リップルタンク内の波の高さの測定など
素粒子物理：画像検出器に入射する粒子の分布の測定など
2次元データは、これまでのグラフで表示されてきたような1次元のリストよりも視覚化が難しいですが、そのような場合には密度プロットが役に立ちます。これは、色または明るさを使用して値を表す2次元プロットです。
moodleから circular.txt をダウンロードしてjupter notebookのディレクトリにおいてください。
[2]:
import matplotlib.pyplot as plt
import numpy as np
data = np.loadtxt("circular.txt",float)
print(np.shape(data))
(501, 501)
読み込んだデータは、501 × 501 個の２次元データとなっています。
読み込んだデータを２次元プロットするには imshow を使います。
[3]:
plt.imshow(data)
plt.show()
縦軸と横軸の数字はリストのインデックスになっています。またデフォルトでは原点が左上にあります。
原点の位置を左下に設定します。
[4]:
plt.imshow(data, origin="lower")
plt.show()
このままでは色の意味がわからないので カラーバー をつけます。
[5]:
plt.imshow(data, origin="lower")
plt.colorbar()
plt.show()
グレースケールで図示するには以下のようにします。
[6]:
plt.imshow(data, origin="lower")
plt.colorbar()
plt.gray()
plt.show()
計算結果を２次元プロットする¶
上で使用した circular.txt は、\(\cos(\sqrt{x^2 + y^2})\) の計算結果です。次の例では、実際に計算を行ったデータについて二次元プロットを行います。
必要なライブラリをインポートします。
[7]:
import matplotlib.pyplot as plt
from math import sqrt,sin,cos,pi
from numpy import empty
波の高さデータを格納する配列(501 × 501)をつくります。
[8]:
xi = empty([501,501],float)
配列に格納するデータを計算します。
[9]:
for i in range(501):
y = 0.02*i
for j in range(501):
x = 0.02*j
r = sqrt((x)**2+(y)**2)
xi[i,j] = cos(r)
２次元プロットを作成します。
[12]:
plt.imshow(xi, origin="lower")
plt.colorbar()
plt.viridis() #デフォルトのカラーマップ
plt.show()
[14]:
plt.imshow(xi, origin="lower", extent=(0,10,0,10))
plt.colorbar()
plt.viridis() #デフォルトのカラーマップ
plt.show()
課題１：波の干渉¶
池に小石を落とすと、落ちた場所から波が放射状に広がります。一定の時間が経った後の波の高さは、一様な円形に広がる正弦波を使って表現することができます。
円の中心が \(x_1\), \(y_1\) にある場合、点 \(x, y\) から円の中心までの距離 \(r_1\)は、
正弦波の高さは、
ここで\(h_0\) は波の振幅、\(k\) は波動ベクトルで、波長 \(\lambda\)を用いて\(k = 2\pi/\lambda\) と表されます。
ここで池にもう一つの小石を落として、別の円形を作成します。
波長と振幅が同じで、中心が\(x_2\), \(y_2\)にある第二の波は、
となります。
次に、波を線形に加算します（これは、波が大きすぎないことを条件としています。）。
\(x, y\)における合成派は、
池の領域を 1 m 四方とし、波の高さの画像を作成するプログラムを作成してください。このとき波長は \(\lambda\) = 5 cm、振幅 1 cm、二つの円の中心は 20 cm 離して設定してください。
イメージを作るには、以下の点に注意してください。
各格子点の波の高さ h の配列を作成します。
その配列を使用して密度プロットを作成します。
例えば501 × 501 ポイントのグリッドを使用して、1 m の正方形をカバーします。
答えは以下のようになります。
課題２：波の干渉の時間変化¶
波の伝わる速度を\(v\)とすると、中心が\(x_1\), \(y_1\) にある時間\(t\)の正弦波の高さは、
波長と振幅が同じで、中心が\(x_2\), \(y_2\)にある第二の波のt秒後の高さは、
となります。
\(t\) = 0 ~ 10 s の場合の波の高さの画像を作成するプログラムを作成してください。(\(v\)=10 cm/s)
[ ]: