関数をグラフ化

次の関数をグラフ化します。

$$
[
f(x) = \frac{\sin(x)}{x}
]
$$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

import numpy as np
import matplotlib.pyplot as plt

# 関数を定義
def f(x):
    return np.sin(x) / x

# xの範囲を設定
x = np.linspace(-10, 10, 1000)
x = x[x != 0]  # ゼロを除外（分母がゼロになるのを防ぐため）

# 関数をグラフにプロット
plt.figure(figsize=(8, 6))
plt.plot(x, f(x), label='$\\frac{\\sin(x)}{x}$')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Graph of $f(x) = \\frac{\\sin(x)}{x}$')
plt.legend()
plt.grid(True)
plt.show()

これにより、 $ (f(x) = \frac{\sin(x)}{x}) $のグラフが表示されます。

[実行結果]

ソースコード解説

このコードは、関数$ (f(x) = \frac{\sin(x)}{x}) $を定義し、それをグラフ化しています。

1. numpy ライブラリを使って、連続的な$ (x) $の範囲を設定しています。
   linspace 関数は、-10から10までの間を等間隔に分割し、1000点の$ (x) $座標を生成しています。

2. ただし、分母がゼロになる$ (x=0) $の点を除外するために、x から$ (x=0) $の点を取り除いています。

3. f(x) の値を計算し、関数の値を取得します。

4. matplotlib を使ってグラフを描画します。
   生成した$ (x) $座標とそれに対応する関数の値をプロットしています。

5. グラフの軸ラベル、タイトル、凡例、グリッド線などが設定されています。

このコードを実行すると、関数$ (f(x) = \frac{\sin(x)}{x}) $のグラフが生成され、その特性が視覚化されます。

特に、$ (x) $がゼロに近いところでの挙動や、$ (x) $が大きくなるにつれて関数がゼロに収束する様子などが観察できます。

グラフ解説

[実行結果]

このグラフは関数 $ (f(x) = \frac{\sin(x)}{x}) $を表しています。
この関数は、$ (x) $の値がゼロに近づくにつれて、分母がゼロになりますが、 $ (\sin(x)) $は$ (x) $がゼロの近くでもゼロにならないため、関数は$ (x=0) $で定義されます。

グラフは、$ (x) $がゼロに近いところで$ (\frac{\sin(x)}{x}) $の値が大きく振動し、$ (x) $の絶対値が大きくなるにつれて、関数の値がゼロに近づく様子を示しています。

関数は$ (x=0) $で未定義ですが、その周辺で連続に振る舞い、$ (x) $が大きくなると関数の値はゼロに収束していきます。

グラフ全体が周期的な特徴を持ち、$ (x) $がゼロに近いところで特に振動が大きくなっています。

生物学的ネットワークの解析

NetworkXはPythonのライブラリで、グラフ理論を操作するのに便利です。

生物学的な問題を解決するために、例えば、タンパク質間の相互作用ネットワークや遺伝子の相互作用ネットワークなど、生物学的ネットワークの解析に使用することができます。

以下は、NetworkXを使用してランダムな生物学的ネットワークを作成し、そのネットワークを可視化する例です。

このコードでは、ネットワークのノードはタンパク質や遺伝子を表し、エッジはそれらの間の相互作用を示します。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

import networkx as nx
import matplotlib.pyplot as plt
import random

# ランダムな生物学的ネットワークを生成するための仮のデータ生成
num_nodes = 20  # ノード数
num_edges = 30  # エッジ数

# 空の有向グラフを作成
G = nx.DiGraph()

# ランダムなノードを追加
for i in range(num_nodes):
    G.add_node(f'Node_{i}')

# ランダムなエッジを追加
for _ in range(num_edges):
    node1 = random.choice(list(G.nodes))
    node2 = random.choice(list(G.nodes))
    if node1 != node2:  # 自己ループを避けるため
        G.add_edge(node1, node2)

# グラフの描画
plt.figure(figsize=(8, 6))
pos = nx.spring_layout(G)  # レイアウトの設定
nx.draw(G, pos, with_labels=True, node_size=300, node_color='skyblue', font_weight='bold', font_size=8, edge_color='gray')
plt.title('Random Biological Network')
plt.show()

この例では、ランダムなノードとエッジを持つ生物学的ネットワークを作成して描画しています。

実際の生物学的ネットワークの場合、ネットワーク構造やデータの特性に合わせて異なる手法やアプローチが必要になるかもしれません。

[実行結果]

ソースコード解説

このコードは、NetworkXとMatplotlibを使用してランダムな生物学的ネットワークを生成し可視化するものです。

1. ライブラリのインポート:

networkxとmatplotlib.pyplotを使っています。
randomモジュールも使われています。

2. ネットワークの生成:

num_nodesで指定された数のノードとnum_edgesで指定された数のエッジを持つ有向グラフを作成します。
まず、空の有向グラフGを作成し、20個のノードをNode_0からNode_19まで追加します。

3. エッジの追加:

num_edgesの数だけランダムなエッジを追加します。
ランダムに2つのノードを選び、それらを結ぶエッジを有向グラフに追加します。
ただし、選ばれた2つのノードが同じでないことを確認して自己ループを避けます。

4. グラフの可視化:

plt.figure()で新しい図を作成し、nx.spring_layout()でノードのレイアウトを設定します。
そして、nx.draw()を使用してノードとエッジを可視化します。
ノードはskyblueの色で、ラベルも表示されます。
エッジはgrayの色で描かれます。
最後にplt.show()でグラフが表示されます。

このコードは、NetworkXを使ってランダムな生物学的ネットワークを生成し、その構造を視覚的に表現しています。

これにより、ネットワーク内のノードとエッジのつながりを視覚化し、ネットワーク全体の特性を把握することができます。

結果解説

[実行結果]

このコードは、NetworkXを使用して作成されたランダムな生物学的ネットワークを可視化しています。
ネットワークは20個のノードと30個のエッジで構成されています。

ノード（Node）:

各ノードはタンパク質や遺伝子を表しており、”Node_0”、”Node_1”などの名前で識別されています。
この例では、ネットワーク内のノードは20個生成されています。

エッジ（Edge）:

ノード間の相互作用を示しています。
ここでは30個のランダムなエッジがあり、ノード間の結びつきを表しています。
エッジは有向グラフとして表現されており、ノード間の特定の方向性を持っています。

グラフの表示:

nx.spring_layout()を使用してグラフのレイアウトを設定し、ノードの配置を行っています。
ノードのラベルや色、エッジの色などが設定されています。
グラフ全体がランダムな配置で描かれ、ノード同士のつながりが視覚的に表現されています。

このようなネットワークの可視化は、ノード間の関係性やネットワーク全体の特性を理解し、分析するための手助けとなります。

生物学的なネットワークの場合、実際のデータを解析し、そのネットワークの特徴や機能を理解するのに役立ちます。

シュレーディンガー方程式

シュレーディンガー方程式は、量子力学における基本的な方程式の一つであり、波動関数を使って系の時間発展を記述します。

１次元の自由粒子のシュレーディンガー方程式は以下のように表されます。

$$
[ -\frac{\hbar^2}{2m} \frac + V(x)\psi = E\psi ]
$$

ここで、$ (\hbar) $はプランク定数の割合、$ (m) $は粒子の質量、$ (\psi) $は波動関数、$ (V(x)) $はポテンシャル関数、$ (E) $はエネルギーです。

この方程式を解析的に解くことは多くの場合難しいですが、数値的に解くことは可能です。

Pythonで数値解を求め、グラフ化するコードを示します。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh_tridiagonal

# パラメータ設定
h_bar = 1.0
m = 1.0
omega = 1.0  # 角振動数の設定
N = 1000
L = 10.0
x = np.linspace(0, L, N)
dx = x[1] - x[0]

# ポテンシャル関数を定義
def potential_function(x):
    return 0.5 * m * (omega**2) * x**2  # 調和振動子ポテンシャル

# 対角成分と非対角成分を計算
diagonal = np.ones(N) / dx**2 + potential_function(x)
off_diagonal = -0.5 * np.ones(N - 1) / dx**2

# トリディアゴナル行列を作成して固有値を求める
eigenvalues, eigenvectors = eigh_tridiagonal(diagonal, off_diagonal)

# 波動関数のグラフ化
plt.figure(figsize=(8, 6))
for i in range(4):  # 最初の４つの基底状態の波動関数を描画
    psi = eigenvectors[:, i]
    plt.plot(x, psi, label=f'n = {i+1}')

plt.title('Wavefunctions of the Quantum Harmonic Oscillator')
plt.xlabel('x')
plt.ylabel('psi(x)')
plt.legend()
plt.grid(True)
plt.show()

このコードでは、単純な調和振動子ポテンシャルを持つシュレーディンガー方程式の基底状態の波動関数を計算し、プロットしています。

これにより、最初の４つの基底状態の波動関数が描画されます。

[実行結果]

ソースコード解説

このコードは、量子調和振動子のポテンシャルに対するシュレーディンガー方程式の固有値問題を解き、最初の４つの基底状態の波動関数をプロットしています。

1. パラメータ設定:

• $ ( \hbar ) $: プランク定数
• $ ( m ) $: 粒子の質量
• $ ( \omega ) $: 角振動数
• $ ( N ) $: 空間分割の数
• $ ( L ) $: 空間の長さ
• $ ( x ) $: 空間の区間

2. ポテンシャル関数:

potential_function(x)で定義され、量子調和振動子ポテンシャル$ ( \frac{1}{2} m \omega^2 x^2 ) $を表します。

3. 対角成分と非対角成分の計算:

トリディアゴナル行列の対角成分と非対角成分を計算します。

4. eigh_tridiagonal関数:

SciPyのeigh_tridiagonal関数を使用して、トリディアゴナル行列の固有値と固有ベクトルを求めます。
これらは波動関数と対応するエネルギー固有値です。

5. 波動関数のプロット:

最初の４つの基底状態の波動関数をforループで計算し、matplotlibを使用してそれらをグラフ化しています。
それぞれの波動関数は n = 1, 2, 3, 4 とラベル付けされており、量子調和振動子の波動関数の性質を示しています。

結果解説

[実行結果]

このグラフは、調和振動子ポテンシャルを持つシュレーディンガー方程式の基底状態の波動関数を表しています。

調和振動子ポテンシャルは、$ (\frac{1}{2}m\omega^2x^2) $と表され、量子力学で非常に重要なポテンシャルです。

プロットされた曲線は、それぞれの波動関数を表しています。
グラフの x 軸は位置を示し、y 軸は各波動関数の値を表しています。
波動関数は、量子力学における粒子の振る舞いを表現し、特定のエネルギー状態で粒子がどのように分布するかを示します。

この例では、最初の４つの基底状態の波動関数がプロットされています。
これらの波動関数は、それぞれ異なるエネルギー状態を表しており、基底状態から順に高次のエネルギー状態が表示されています。

波動関数はゼロ点振動（原点で振幅がゼロ）から始まり、対称性を持って左右に広がります。
それぞれの波動関数は異なる空間的な特性を持っており、基底状態から高次のエネルギー状態へとエネルギーが増加するにつれて、振動の幅が大きくなります。

ローレンツ方程式

非線形微分方程式であるローレンツ方程式を解き、その結果を3次元のグラフで表示してみましょう。

ローレンツ系は次の3つの非線形連立微分方程式からなります。

$$
dx/dt = σ(y - x)
$$
$$
dy/dt = x(ρ - z) - y
$$
$$
dz/dt = xy - βz
$$

ここで、$σ、ρ、β$は正の定数です。

以下のPythonコードは、上記のローレンツ方程式を解き、その結果を3Dグラフで表示します。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def lorenz_system(current_state, t):
    sigma = 10.0 # σの値
    rho = 28.0 # ρの値
    beta = 8.0 / 3.0 # βの値

    x, y, z = current_state # 現在のx, y, zの値
    dx_dt = sigma * (y - x) # dx/dt
    dy_dt = x * (rho - z) - y # dy/dt
    dz_dt = x * y - beta * z # dz/dt

    return [dx_dt, dy_dt, dz_dt]

# 初期状態
initial_state = [1.0, 1.0, 1.0]

# 時間の範囲を定義
t = np.arange(0.0, 50.0, 0.01)

# 微分方程式を解く
solution = odeint(lorenz_system, initial_state, t)

# 結果を3Dグラフで表示
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(solution[:, 0], solution[:, 1], solution[:, 2])
ax.set_title('Lorenz Attractor')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()

このコードを走らせると、解の軌道がローレンツアトラクターと呼ばれる蝶形の図形を形成します。

これはカオス理論を象徴する図形としてよく知られています。

[実行結果]

ソースコード解説

このコードは、ローレンツ方程式を数値的に解き、結果を3Dグラフとして表示するものです。

1. 必要なライブラリをインポートする:

numpy（数値計算のため）、scipy.integrateからのodeint（常微分方程式を解くため）、matplotlib.pyplot（グラフ描画のため）、mpl_toolkits.mplot3dからのAxes3D（3Dプロットのため）をインポートします。

2. ローレンツ方程式を定義する:

lorenz_system関数は、ローレンツ方程式を定義しています。
この方程式は3つの微分方程式からなり、特定のパラメータ（sigma、rho、beta）と現在の状態（x、y、z）を受け取ります。

それぞれの微分方程式は次のように表されます。
$$
dx/dt = sigma * (y - x)
$$
$$
dy/dt = x * (rho - z) - y
$$
$$
dz/dt = x * y - beta * z
$$

3. 初期状態を設定する:

initial_stateは、$x、y、z$の初期値を$1.0$として設定します。

4. 時間の範囲を定義する:

tは$0$から$50$までの時間の範囲を$0.01$刻みで定義します。

5. 微分方程式を解く:

odeint関数を使用して、lorenz_system関数に初期状態と時間の範囲を渡して微分方程式を解きます。

6. 結果を3Dグラフで表示する:

solutionは微分方程式の解です。
これを3Dグラフでプロットして、X、Y、Z軸に解の値をプロットし、Lorenz Attractorというタイトルをつけます。

これにより、ローレンツアトラクタとして知られる複雑な非線形系の挙動が視覚化されます。

3次元ガウス関数は、多変量ガウス分布を表現するための関数です。

以下は3次元ガウス関数の式です。

$$
[ f(x, y) = \frac{1}{2\pi\sigma_x\sigma_y}e^{-\frac{1}{2}\left(\frac{(x-\mu_x)^2}{\sigma_x^2} + \frac{(y-\mu_y)^2}{\sigma_y^2}\right)} ]
$$

ここで、$ (x) $と$ (y) $は変数で、$ (\mu_x) $と$ (\mu_y) $はそれぞれ$ (x) $と$ (y) $の平均、$ (\sigma_x) $と$ (\sigma_y) $は標準偏差です。

Pythonで3次元ガウス関数をグラフ化するための例を示します。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D

def gaussian_2d(x, y, mux, muy, sigmax, sigmay):
    return (1. / (2 * np.pi * sigmax * sigmay) * 
            np.exp(-0.5 * ((x - mux)**2 / sigmax**2 + (y - muy)**2 / sigmay**2)))

# ガウス関数のパラメータ設定
mux = 0
muy = 0
sigmax = 1
sigmay = 1

x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
x, y = np.meshgrid(x, y)
z = gaussian_2d(x, y, mux, muy, sigmax, sigmay)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x, y, z, cmap='viridis')

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('3D Gaussian Function')

plt.show()

このコードでは、2つの変数$ (x) $と$ (y) $のガウス関数を計算し、3次元プロットを作成しています。

平均と標準偏差の値を変更すると、関数の形状がどのように変わるかを観察できます。

[実行結果]

ソースコード解説

このコードは、2変数の2次元ガウス関数を作成し、その関数を3次元プロットで視覚化するためのものです。

以下、詳細な説明です：

1. import文:

• matplotlib.pyplotからpltをインポートしています。
  これはグラフの描画に使用されます。
• numpyからnpをインポートしています。
  これは数値計算に使用されます。
• mpl_toolkits.mplot3dからAxes3Dをインポートしています。
  これは3次元プロットのための機能を提供します。

2. gaussian_2d関数:

• 2変数の2次元ガウス関数を定義しています。
  この関数は、座標 (x, y) に対するガウス関数の値を返します。

3. ガウス関数のパラメータ設定:

• mux、muy、sigmax、sigmay は、ガウス関数の平均と標準偏差を設定するためのパラメータです。

4. np.linspace()を使って x と y の値を生成:

• np.linspace()は指定された範囲内で一定間隔の値を生成します。

5. np.meshgrid()を使って x と y の格子を作成:

• np.meshgrid()は x と y の2つの1次元配列を受け取り、それぞれの組み合わせに対する格子を作成します。

6. gaussian_2d()関数を使って x と y に対するガウス関数の値 z を計算:

• gaussian_2d()関数は、先ほど定義したガウス関数を使用して、各 (x, y) 座標でのガウス関数の値を計算します。

7. plt.figure()を使って図を作成:

• plt.figure()は新しい図を作成します。

8. fig.add_subplot()を使って3Dサブプロットを追加:

• fig.add_subplot()は図に新しい3次元のサブプロットを追加します。

9. ax.plot_surface()を使って3Dプロットを作成:

• ax.plot_surface()は x、y、z 座標上のデータから3次元プロットを作成します。
  cmap='viridis'はカラーマップを設定します。

10. ax.set_xlabel()、ax.set_ylabel()、ax.set_zlabel()を使って軸ラベルを設定:

• それぞれ X、Y、Z 軸のラベルを設定します。

11. ax.set_title()を使ってグラフタイトルを設定:

• グラフにタイトルを追加します。

12. plt.show()でグラフを表示:

• plt.show()は作成したグラフを表示します。

このコードは、2変数の2次元ガウス関数を計算し、それを3次元プロットとして視覚化しています。

球体の方程式

3次元方程式の例として、球体の方程式を挙げて、そのグラフ化を行います。

球体の方程式は次のように表されます。

$$
[ x^2 + y^2 + z^2 = r^2 ]
$$

これは、原点を中心とし、半径が$ (r) $の球の方程式です。

Pythonでこの方程式をグラフ化するには、Matplotlibを使用して3次元のプロットを行います。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

r = 1  # 球の半径

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# x, y, zの値を生成
theta = np.linspace(0, 2 * np.pi, 100)
phi = np.linspace(0, np.pi, 100)
x = r * np.outer(np.cos(theta), np.sin(phi))
y = r * np.outer(np.sin(theta), np.sin(phi))
z = r * np.outer(np.ones(100), np.cos(phi))

# 球をプロット
ax.plot_surface(x, y, z, color='b', alpha=0.5)

# 軸ラベルの設定
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')

plt.show()

このコードでは、球の方程式をパラメータ化して球体を描画しています。

球の半径や解像度を変更することで、異なる球を描画することができます。

[実行結果]

ソースコード解説

このコードは、Matplotlibを使用して3次元の球を描画するものです。

それぞれの部分について説明します。

1. モジュールのインポート:

1
2
3

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

必要なライブラリをインポートしています。
matplotlib.pyplotはグラフを描画するために使用され、mpl_toolkits.mplot3dは3次元プロットをサポートするために必要です。
numpyは数学的な演算を行うために使用されます。

2. 球の半径の設定:

1

r = 1  # 球の半径

描画する球の半径を設定しています。

3. FigureとAxesの設定:

1
2

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

plt.figure()を使用してFigureオブジェクトを作成し、fig.add_subplot()でAxesオブジェクトを作成しています。
ここでは、3次元プロットを行うためにprojection='3d'を指定しています。

4. x, y, zの値の生成:

1
2
3
4
5

theta = np.linspace(0, 2 * np.pi, 100)
phi = np.linspace(0, np.pi, 100)
x = r * np.outer(np.cos(theta), np.sin(phi))
y = r * np.outer(np.sin(theta), np.sin(phi))
z = r * np.outer(np.ones(100), np.cos(phi))

球の表面を構成する$x, y, z$座標の値を生成しています。
np.linspace()は指定された範囲で等間隔の値を生成し、np.outer()は2つの配列の外積を計算します。

5. 球のプロット:

1

ax.plot_surface(x, y, z, color='b', alpha=0.5)

ax.plot_surface()を使用して球を描画しています。
ここでは、生成したx, y, z座標の値を使用し、color='b'で青色、alpha=0.5で透明度を設定しています。

6. 軸ラベルの設定:

1
2
3

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')

$x$軸、$y$軸、$z$軸のラベルを設定しています。

7. グラフの表示:

1

plt.show()

最後にplt.show()を使ってグラフを表示します。
これにより、3次元の球が表示されます。

ロジスティックマップ

非線形の関数としてロジスティックマップを例に挙げます。

ロジスティックマップは、次のような非線形差分方程式で定義されます：

$$
x[n+1] = r * x[n] * (1 - x[n])
$$

ここで、$ n $は離散時間、$ x[n] $は時刻nにおける値、$ r $は制御パラメータです。

この方程式は生物個体数などの成長モデルに用いられることがあります。

パラメータ$r$により、周期$ 2, 4, 8, …, caos $など様々なダイナミクスを見ることができます。

以下に、ロジスティックマップをPythonでグラフ化するコードを示します。

これは$r$を$0.5$から$4.0$まで$0.01$間隔で変化させ、$x$の時間発展をプロットするものです：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

import numpy as np
import matplotlib.pyplot as plt

def logistic(r, x):
    return r * x * (1 - x)

x = np.linspace(0, 1)
fig, ax = plt.subplots(1, 1, figsize=(6, 10), dpi=80)
for r in np.linspace(0.5, 4.0, 400):
    X = np.empty(1000)
    X[0] = 0.1
    for i in range(999):
        X[i+1] = logistic(r, X[i])
    ax.plot([r]*1000, X, ',k', alpha=0.2)

ax.set_xlim(0.4, 4)
ax.set_title("Bifurcation diagram")
plt.show()

これによりビフルケーション図が得られます。

点々とした部分がカオス性を示し、この隙間の生じるパターンはフラクタル性を示します。

また、一定範囲の$ r $での値の振動はピリオドダブリングと呼ばれ、カオスへと向かう典型的な挙動を表しています。

ソースコード解説

このPythonのコードは、ロジスティック写像の分岐図（Bifurcation diagram）を描画するものです。

それぞれの部分について詳しく説明します。

ライブラリのインポート

1
2

import numpy as np
import matplotlib.pyplot as plt

• numpyは数値計算をサポートするPythonライブラリで、配列操作や数値計算を高速に行います。
• matplotlib.pyplotはグラフ描画ライブラリです。

ロジスティック写像の関数

1
2

def logistic(r, x):
    return r * x * (1 - x)

• logistic関数はロジスティック写像の式を表しています。

分岐図の作成

1
2
3
4
5
6
7
8

x = np.linspace(0, 1)
fig, ax = plt.subplots(1, 1, figsize=(6, 10), dpi=80)
for r in np.linspace(0.5, 4.0, 400):
    X = np.empty(1000)
    X[0] = 0.1
    for i in range(999):
        X[i+1] = logistic(r, X[i])
    ax.plot([r]*1000, X, ',k', alpha=0.2)

• np.linspace(0, 1)で0から1までの間を等間隔に分割した配列xを作成します。
• fig, ax = plt.subplots(1, 1, figsize=(6, 10), dpi=80)でグラフのサイズを設定します。
• for r in np.linspace(0.5, 4.0, 400):でパラメータrの範囲を指定します。
  rは0.5から4.0の間を400ステップで変化します。
• ロジスティック写像を使って、与えられたパラメータrにおける1000ステップの反復を計算し、それをプロットします。
• ax.plot([r]*1000, X, ',k', alpha=0.2)で分岐図を描画します。
  [r]*1000でrを1000回繰り返した配列を作成し、それとXをプロットします。

グラフの設定

1
2
3

ax.set_xlim(0.4, 4)
ax.set_title("Bifurcation diagram")
plt.show()

• ax.set_xlim(0.4, 4)でx軸の範囲を設定します。
• ax.set_title("Bifurcation diagram")でグラフのタイトルを設定します。
• plt.show()でグラフを表示します。

このコードは、ロジスティック写像による分岐図を描画するためのもので、特定のパラメータ範囲での振る舞いを可視化します。

流体力学

流体力学の問題の1つとして、流れの中の速度分布を考えてみましょう。

例として、流れが円管内を通る場合の速度分布を解析してみます。

ナビエ・ストークス方程式を仮定し、理想的な流れを対象にします。

以下は、円管内の流れにおける速度分布の解析とグラフ化の例です。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

# パラメータの設定
R = 1.0  # 円管の半径
viscosity = 0.1  # 粘度係数

# ナビエ・ストークス方程式の関数
def equation(u):
    return u - (1 / R) * (0.5 * viscosity) * ((1 - (u ** 2)) / u)

# 初期推定値
initial_guess = 0.1

# 方程式を解いて速度分布を求める
u_solution = fsolve(equation, initial_guess)

# 半径方向の位置
r_values = np.linspace(0, R, 100)

# 速度分布の計算
u_values = (1 / R) * (0.5 * viscosity) * ((1 - (u_solution ** 2)) / u_solution) * (1 - (r_values / R) ** 2)

# 結果をグラフ化
plt.figure(figsize=(8, 6))
plt.plot(r_values, u_values)
plt.xlabel('Radius')
plt.ylabel('Velocity')
plt.title('Velocity Distribution in Pipe Flow')
plt.grid(True)
plt.show()

このコードでは、ナビエ・ストークス方程式の近似解析を行い、円管内の速度分布を計算しています。

結果はMatplotlibを使用してグラフ化されます。

得られた速度分布は、円管内の中心から周辺への速度変化を示しています。

[実行結果]

ソースコード解説

このPythonのコードは、円管内の流れにおける速度分布を解析し、結果をグラフ化するためのものです。

それぞれの部分について詳しく説明します。

ライブラリのインポート

1
2
3

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

• numpyは数値計算をサポートするPythonライブラリで、行列や配列などの高性能な数値計算機能を提供します。
• matplotlib.pyplotはデータを視覚化するためのグラフ描画ライブラリです。
  pltとして一般的にエイリアスされます。
• scipy.optimize.fsolveは非線形方程式の数値解を見つけるためのSciPyライブラリの関数です。

パラメータの設定

1
2

R = 1.0  # 円管の半径
viscosity = 0.1  # 粘度係数

• Rは円管の半径を表し、この例では$1.0$に設定されています。
• viscosityは流体の粘度係数を表し、この例では$0.1$に設定されています。

ナビエ・ストークス方程式の関数

1
2

def equation(u):
    return u - (1 / R) * (0.5 * viscosity) * ((1 - (u ** 2)) / u)

• equation関数は、円管内の流れにおける速度分布を表すナビエ・ストークス方程式の非線形方程式を定義します。

方程式を解いて速度分布を求める

1
2

initial_guess = 0.1
u_solution = fsolve(equation, initial_guess)

• fsolve関数を使用して、equation関数を解くことで速度分布の解析解を計算します。
  initial_guessは初期推定値を表します。

半径方向の位置の計算

1

r_values = np.linspace(0, R, 100)

• np.linspace関数を使用して、$0$から円管の半径（R）までの範囲を$100$等分した半径方向の位置を作成します。

速度分布の計算

1

u_values = (1 / R) * (0.5 * viscosity) * ((1 - (u_solution ** 2)) / u_solution) * (1 - (r_values / R) ** 2)

• ナビエ・ストークス方程式の解から、円管内の各位置での速度分布を計算します。

結果をグラフ化

1
2
3
4
5
6
7

plt.figure(figsize=(8, 6))
plt.plot(r_values, u_values)
plt.xlabel('Radius')
plt.ylabel('Velocity')
plt.title('Velocity Distribution in Pipe Flow')
plt.grid(True)
plt.show()

• matplotlibを使用して、計算された速度分布を円管の半径に対する速度としてグラフ化します。
  plt.xlabel、plt.ylabel、plt.titleを使用して、軸のラベルとグラフのタイトルを設定し、plt.grid(True)でグリッドを表示します。
  plt.show()でグラフを表示します。

このコードは、ナビエ・ストークス方程式を使用して円管内の流れにおける速度分布を解析し、その結果を可視化するものです。

グラフ解説

このグラフは円管内の流れにおける速度分布を示しています。

[実行結果]

x軸は円管の半径を示し、y軸はその位置における流速を表しています。

解析した速度分布は、理想的な流れ条件を仮定しています。
円管内の中心$ （r=0）$では、流速は最大になります。
そして、円管の壁に近づくにつれて流速はゼロに近づきます。

この結果は、円管内の流れにおいて、中心部分の流れが最も速く、壁に近づくほど流速が低くなることを示しています。
このような速度分布は、レイノルズ数や流体の粘度などの条件によって変化する可能性がありますが、円管流れの一般的な傾向を示しています。

また、このグラフは円管内の速度分布が放射状に変化していることを示しています。
中心部から外側に向かって流速が減少していることが分かります。
このような速度プロファイルは、流体が壁面で摩擦を受けることによって生じる、粘性効果によるものです。

グループ化／集約

Polarsの一つの強みは、大規模データセットを高速に処理する能力です。

その最適化された処理エンジンのおかげで、大量のデータを即座に変換、整理したり、集計することが容易となり、パフォーマンスを大いに向上させます。

また、シンプルで明確なAPIもPolarsの大きな魅力の一つです。

これにより、ユーザーはpandasのような一般的なデータ分析ライブラリと同様のユーザビリティを享受しながら、より大きなデータ量を扱うことができます。

以下に、Polarsを使ったデータのグループ化と集約のサンプルコードを示します。

1
2
3
4
5
6
7
8
9
10
11
12
13

import polars as pl

# データフレームの作成
df = pl.DataFrame({
    "fruits": ["apple", "banana", "banana", "apple", "apple", "banana", "banana", "apple"],
    "sales": [5, 6, 8, 9, 5, 6, 7, 9],
    "city": ["Tokyo", "Tokyo", "Tokyo", "London", "London", "London", "New York", "New York"]
})

# 都市ごと、果物ごとの総売上を計算
df_grouped = df.groupby(["city", "fruits"]).agg([pl.col("sales").sum().alias("total_sales")])

print(df_grouped)

上記のコードは、都市ごとに果物の売上合計を計算するシンプルな例です。

都市と果物ごとにグループ化し、各グループの“sales”列の合計を計算します。

[実行結果]

shape: (6, 3)
┌──────────┬────────┬─────────────┐
│ city     ┆ fruits ┆ total_sales │
│ ---      ┆ ---    ┆ ---         │
│ str      ┆ str    ┆ i64         │
╞══════════╪════════╪═════════════╡
│ Tokyo    ┆ apple  ┆ 5           │
│ Tokyo    ┆ banana ┆ 14          │
│ New York ┆ banana ┆ 7           │
│ London   ┆ apple  ┆ 14          │
│ London   ┆ banana ┆ 6           │
│ New York ┆ apple  ┆ 9           │
└──────────┴────────┴─────────────┘

shape: (6, 3) はこのデータフレームが6行3列であることを示しています。
列は city、fruits、そして total_sales となります。

それぞれの行を見てみましょう：

• 1行目：
  東京でのリンゴの総売上は5と計算されています。
• 2行目：
  東京でのバナナの総売上は14と計算されています。
• 3行目：
  ニューヨークでのバナナの総売上は7と計算されています。
• 4行目：
  ロンドンでのリンゴの総売上は14と計算されています。
• 5行目：
  ロンドンでのバナナの総売上は6と計算されています。
• 6行目：
  ニューヨークでのリンゴの総売上は9と計算されています。

これにより、各都市で各フルーツの総販売量を直観的に理解することができます。

このようなデータの表示は、日々の販売活動の理解や、将来的な販売戦略の策定に非常に有用です。

フルーツの種類と都市に基づいてデータを分析することで、どの商品がどの市場でよく売れているのか洞察を得ることが可能です。

Polarsの概要

Polarsは、高速なDataFrame操作を可能にするライブラリで、Apache Arrowの列指向メモリモデルを使用しています。

Polarsは、Rustで書かれており、多スレッドのクエリエンジンを使用して効率的な並列処理を実現しています。

また、Polarsはオープンソースであり、MITライセンスで利用できます。

サンプルソース

以下に、Polarsを使用してDataFrameを作成し、そのDataFrameに対して操作を行う例を示します。

まず、必要なライブラリをインポートします。

1

import polars as pl

次に、DataFrameを作成します。

1
2
3
4
5
6
7

df = pl.DataFrame(
 {
     "name": ["John", "Sue", "Tom", "Alice"],
     "age": [23, 21, 19, 18],
     "grade": [90, 85, 78, 92],
 }
)

次に、DataFrameに対して操作を行います。

例えば、特定の列を選択したり、特定の条件に基づいて行をフィルタリングしたり、新しい列を作成したりすることができます。

1
2
3
4
5
6
7
8

# 特定の列を選択
selected_df = df.select(["name", "age"])

# 特定の条件に基づいて行をフィルタリング
filtered_df = df.filter(pl.col("age") > 20)

# 新しい列を作成
df = df.with_columns(pl.col("age") + 1)

最後に、DataFrameを表示します。

1
2
3

print(selected_df)
print(filtered_df)
print(df)

[実行結果]

shape: (4, 2)
┌───────┬─────┐
│ name  ┆ age │
│ ---   ┆ --- │
│ str   ┆ i64 │
╞═══════╪═════╡
│ John  ┆ 23  │
│ Sue   ┆ 21  │
│ Tom   ┆ 19  │
│ Alice ┆ 18  │
└───────┴─────┘
shape: (2, 3)
┌──────┬─────┬───────┐
│ name ┆ age ┆ grade │
│ ---  ┆ --- ┆ ---   │
│ str  ┆ i64 ┆ i64   │
╞══════╪═════╪═══════╡
│ John ┆ 23  ┆ 90    │
│ Sue  ┆ 21  ┆ 85    │
└──────┴─────┴───────┘
shape: (4, 3)
┌───────┬─────┬───────┐
│ name  ┆ age ┆ grade │
│ ---   ┆ --- ┆ ---   │
│ str   ┆ i64 ┆ i64   │
╞═══════╪═════╪═══════╡
│ John  ┆ 24  ┆ 90    │
│ Sue   ┆ 22  ┆ 85    │
│ Tom   ┆ 20  ┆ 78    │
│ Alice ┆ 19  ┆ 92    │
└───────┴─────┴───────┘

このコードは、Polarsを使用してDataFrameを作成し、そのDataFrameに対して操作を行い、その結果を表示します。

このように、Polarsを使用することで、DataFrameの操作を効率的に行うことができます。

Polarsは、高速なDataFrame操作を可能にするライブラリで、多くの操作が可能です。

以下に、一部の操作を示します。

1. 行の選択：

DataFrame内の1行を選択するには、row()メソッドを使用して行番号を渡します。
複数の行を選択するためには、filter()関数を使用します。

1
2

# 最初の行を取得する
df.row(0)

[実行結果]

('John', 24, 90)

1
2
3

# 'name'列の値が'Tom'の行を取得する
res = df.filter(pl.col('name') == 'Tom')
print(res)

[実行結果]

shape: (1, 3)
┌──────┬─────┬───────┐
│ name ┆ age ┆ grade │
│ ---  ┆ --- ┆ ---   │
│ str  ┆ i64 ┆ i64   │
╞══════╪═════╪═══════╡
│ Tom  ┆ 20  ┆ 78    │
└──────┴─────┴───────┘

2. 論理演算子の使用：

Polarsでは、論理演算子を使用して複数の条件を指定できます。

1
2
3

# 'name'列の値が'John'または'Alice'の行を取得する
res = df.filter((pl.col('name') == 'John') | (pl.col('name') == 'Alice'))
print(res)

[実行結果]

shape: (2, 3)
┌───────┬─────┬───────┐
│ name  ┆ age ┆ grade │
│ ---   ┆ --- ┆ ---   │
│ str   ┆ i64 ┆ i64   │
╞═══════╪═════╪═══════╡
│ John  ┆ 24  ┆ 90    │
│ Alice ┆ 19  ┆ 92    │
└───────┴─────┴───────┘

3. 結果の格納：

Polarsは演算結果を新しい列として格納する際の記述方法に優れています。
例えば、新しい列を作成するときには、普通の四則演算以外の結果も格納できます。

1
2
3

# 'age'列の値を2で割った結果を新しい列'age_half'に格納する
df2 = df.with_columns(pl.col("age") / 2)
print(df2)

[実行結果]

shape: (4, 3)
┌───────┬──────┬───────┐
│ name  ┆ age  ┆ grade │
│ ---   ┆ ---  ┆ ---   │
│ str   ┆ f64  ┆ i64   │
╞═══════╪══════╪═══════╡
│ John  ┆ 12.0 ┆ 90    │
│ Sue   ┆ 11.0 ┆ 85    │
│ Tom   ┆ 10.0 ┆ 78    │
│ Alice ┆ 9.5  ┆ 92    │
└───────┴──────┴───────┘

これらの操作は、Polarsの強力な機能の一部を示しています。

Polarsは、DataFrameの操作を効率的に行うための強力なツールであり、これらの操作を通じて、データ分析やデータ操作の幅広いタスクを効率的に実行することができます。