ランベルト反射（Lambertian reflection）

乱反射に関する方程式の例として、ランベルト反射（Lambertian reflection）があります。

この反射モデルでは、光が表面に当たった後、ランダムな方向に均等に散乱されると仮定されます。

このモデルでは、反射率は入射角に依存せず、すべての方向に対して均一です。

ランベルト反射の方程式は以下のように表されます:

$$
L_o(\omega_o) = \frac{c}{\pi} L_i(\omega_i) \max(0, \omega_i \cdot n)
$$

ここで、$ ( L_o(\omega_o) ) $は出射光の放射束密度、$ ( L_i(\omega_i) ) $は入射光の放射束密度、$ ( \omega_o ) $は視線の方向、$ ( \omega_i ) $は入射光の方向、$ ( n ) $は法線ベクトル、$ ( c ) $は反射率です。

Pythonでこの式を解き、グラフ化するためのサンプルコードを以下に示します。

このコードでは、ランダムな入射光の方向を仮定し、その方向に応じてランベルト反射による出射光の放射束密度を計算し、その結果をグラフ化します。

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

import numpy as np
import matplotlib.pyplot as plt

def lambertian_reflection(cos_theta_i):
    return np.maximum(0, cos_theta_i) / np.pi

# 入射光の角度を生成
num_samples = 1000
cos_theta_i = np.random.uniform(low=-1, high=1, size=num_samples)
theta_i = np.arccos(cos_theta_i)

# ランベルト反射による放射束密度を計算
lo = lambertian_reflection(cos_theta_i)

# 結果をグラフ化
plt.figure(figsize=(8, 6))
plt.scatter(np.degrees(theta_i), lo, s=5)
plt.title('Lambertian Reflection')
plt.xlabel('Incident Angle (degrees)')
plt.ylabel('Radiant Intensity')
plt.grid(True)
plt.show()

このコードでは、入射光の角度を一様にランダムに生成し、それに対応する出射光の放射束密度を計算しています。

最後に、入射角度に対する出射光の放射束密度をグラフ化しています。

[実行結果]

ソースコード解説

このソースコードは、Pythonを使用してランベルト反射に関する放射束密度を計算し、グラフ化するものです。

以下は、それぞれの部分の詳細な説明です。

インポート:

1
2

import numpy as np
import matplotlib.pyplot as plt

• numpy ライブラリを np としてインポートし、数値計算を行うために使用します。
• matplotlib.pyplot モジュールを plt としてインポートし、グラフを描画するために使用します。

関数定義:

1
2

def lambertian_reflection(cos_theta_i):
    return np.maximum(0, cos_theta_i) / np.pi

• lambertian_reflection 関数は、ランベルト反射による放射束密度を計算する関数です。
• 入射角の余弦を受け取り、ランベルト反射による放射束密度を計算して返します。

入射光の角度生成:

1
2
3

num_samples = 1000
cos_theta_i = np.random.uniform(low=-1, high=1, size=num_samples)
theta_i = np.arccos(cos_theta_i)

• num_samples 個のサンプルを生成します。
• np.random.uniform 関数を使用して、$-1$から$1$までの一様分布から num_samples 個の乱数を生成し、入射光の角度の余弦を生成します。
• np.arccos 関数を使用して、余弦から入射角を計算します。

ランベルト反射による放射束密度の計算:

1

lo = lambertian_reflection(cos_theta_i)

• 先ほど定義した lambertian_reflection 関数を使用して、ランベルト反射による放射束密度を計算します。

結果のグラフ化:

1
2
3
4
5
6
7

plt.figure(figsize=(8, 6))
plt.scatter(np.degrees(theta_i), lo, s=5)
plt.title('Lambertian Reflection')
plt.xlabel('Incident Angle (degrees)')
plt.ylabel('Radiant Intensity')
plt.grid(True)
plt.show()

• グラフの大きさを指定し、plt.figure(figsize=(8, 6)) を使用して新しいグラフを作成します。
• plt.scatter を使用して、散布図をプロットします。
  $x軸$に入射角、$y軸$に放射束密度をプロットします。
• plt.title、plt.xlabel、plt.ylabel を使用して、グラフのタイトルと軸ラベルを設定します。
• plt.grid(True) を使用して、グリッド線を表示します。
• plt.show() を使用して、グラフを表示します。

これらの手順により、ランベルト反射に関する放射束密度の計算結果がグラフ化されます。

結果解説

[実行結果]

上記のPythonコードによって生成されるグラフは、ランベルト反射による放射束密度の振る舞いを示しています。
以下は、グラフの内容の詳細な説明です。

横軸 (Incident Angle):

入射光の角度を示しています。
この角度は度数法で表され、0度が垂直方向、90度が水平方向を意味します。
ランダムに生成された入射光の角度がこの軸に表示されます。

縦軸 (Radiant Intensity):

出射光の放射束密度を示しています。
これはランベルト反射によって得られた出射光の強度を示します。
放射束密度は無次元の値で、ランベルト反射においては最大値が$1$です。

散布図 (Scatter Plot):

各点は、特定の入射角度に対応する出射光の放射束密度を表しています。
この散布図は、入射光の角度に対する出射光の振る舞いを視覚的に示します。
点が高いほど、出射光の強度が大きいことを示します。

タイトル (Title):

グラフのタイトルは「Lambertian Reflection」となっており、このグラフがランベルト反射に関するものであることを示します。

横軸ラベル (X-Axis Label):

「Incident Angle (degrees)」と表示されており、横軸が入射角度を表していることを示します。

縦軸ラベル (Y-Axis Label):

「Radiant Intensity」と表示されており、縦軸が出射光の放射束密度を表していることを示します。

グリッド (Grid):

グラフにはグリッドが表示されており、視覚的にデータの位置を把握しやすくなっています。

このグラフは、入射光の角度に対するランベル反射による出射光の振る舞いを表しており、ランダムな入射光の方向に対してランベル反射がどのように働くかを視覚的に理解するのに役立ちます。

人口増加（ロジスティック成長モデル）

人口増加を表すロジスティック成長モデルを取り上げてみましょう。

このモデルでは、人口が一定の増加率で成長し、限界人口容量に近づくにつれて成長率が減少していくと仮定されます。

ロジスティック成長モデルの方程式は次の通りです：

$$
P(t) = \frac{K}{1 + Ae^{-rt}}
$$

• $ ( P(t) ) $は時間$ ( t ) $における人口数
• $ ( K ) $は限界人口容量
• $ ( A ) $は初期の人口成長率に関連する定数
• $ ( r ) $は成長率の定数

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

import numpy as np
import matplotlib.pyplot as plt

# ロジスティック成長モデルの方程式
def logistic_growth(t, K, A, r):
    return K / (1 + A * np.exp(-r * t))

# パラメータの設定
K = 10000  # 限界人口容量
A = 100    # 初期の人口成長率に関連する定数
r = 0.1    # 成長率の定数

# 時間の範囲を指定
t_values = np.linspace(0, 50, 100)

# 人口数を計算
population = logistic_growth(t_values, K, A, r)

# グラフの描画
plt.figure(figsize=(8, 6))
plt.plot(t_values, population, label='Population')
plt.xlabel('Time')
plt.ylabel('Population')
plt.title('Logistic Growth Model')
plt.grid(True)
plt.legend()
plt.show()

このコードは、ロジスティック成長モデルを解いて人口増加を計算し、時間に対する人口の変化をグラフ化します。

[実行結果]

時間が経つにつれて人口は増加し、限界人口容量に近づくことがグラフからわかります。

ソースコード解説

このソースコードは、Pythonを使用してロジスティック成長モデルに基づいて人口増加をモデル化し、グラフ化するためのものです。

以下ではコードの各部分を詳しく説明します。

1. import numpy as np:

• NumPyライブラリを np としてインポートします。
  NumPyは数値計算を行うためのPythonライブラリであり、配列や行列演算などの機能を提供します。

2. import matplotlib.pyplot as plt:

• Matplotlibの pyplot モジュールを plt としてインポートします。
  Matplotlibはグラフを描画するためのPythonライブラリです。

3. def logistic_growth(t, K, A, r)::

• logistic_growth 関数を定義します。
  この関数は、ロジスティック成長モデルの方程式を実装し、時間$ ( t ) $における人口数を計算します。
• 引数として、時間$ ( t )$、限界人口容量$ ( K )$、初期の人口成長率に関連する定数$ ( A )$、成長率の定数$ ( r ) $を取ります。

4. パラメータの設定:

• K = 10000, A = 100, r = 0.1 のように、ロジスティック成長モデルのパラメータを設定します。
  これらの値はモデルの特性を決定します。

5. t_values = np.linspace(0, 50, 100):

• np.linspace 関数を使って、$0$から$50$までの範囲を$100$の等間隔に分割した配列を生成し、t_values に代入します。
  これは時間の範囲を指定します。

6. population = logistic_growth(t_values, K, A, r):

• logistic_growth 関数を用いて、各時間における人口数を計算し、population に代入します。

7. グラフの描画:

• plt.figure(figsize=(8, 6)) で図のサイズを設定します。
• plt.plot(t_values, population, label='Population') で、時間に対する人口数の曲線をプロットします。
  ラベルは ‘Population’ です。
• plt.xlabel('Time') と plt.ylabel('Population') でx軸とy軸のラベルを設定します。
• plt.title('Logistic Growth Model') でグラフのタイトルを設定します。
• plt.grid(True) でグリッドを表示します。
• plt.legend() で凡例を表示します。
• plt.show() でグラフを表示します。

これにより、ロジスティック成長モデルに基づいて計算された時間と人口数の関係が可視化されます。

結果解説

[実行結果]

このグラフは、ロジスティック成長モデルに基づいて計算された時間と人口数の関係を表しています。

• $x軸$（横軸）は時間を表し、$0$から$50$までの範囲で設定されています。
  時間が経つにつれて、人口の変化を追跡します。
• $y軸$（縦軸）は人口数を表しており、$0$から$10000$の範囲で設定されています。
  これは限界人口容量（K）で、人口がこの値に近づくにつれて成長率が減少します。
• グラフの曲線は、ロジスティック成長モデルに基づいて計算された人口の変化を表しています。
  曲線は時間に対する人口の増加を示し、初めは急速に増加しますが、徐々に成長率が減少し、最終的に限界人口容量に収束します。
• グラフのタイトルは “Logistic Growth Model” となっており、このグラフがロジスティック成長モデルに基づいていることを示しています。
• 凡例には “Population” と表示されており、この曲線が人口数を表していることを示しています。

このグラフからは、時間の経過とともに人口が増加し、限界人口容量に向かって収束する様子が視覚的に理解できます。

ローゼンブロック関数

ローゼンブロック関数は、最適化アルゴリズムの性能テストやベンチマークによく使用される非線形関数です。

ローゼンブロック関数の式は次の通りです：

$$
f(x, y) = (a - x)^2 + b \cdot (y - x^2)^2
$$

ここで、$( a ) と ( b ) $は任意の定数ですが、一般的には$ ( a = 1 ) $および$ ( b = 100 ) $の値が使用されます。

それでは、この関数を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

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

# ローゼンブロック関数の定義
def rosenbrock(x, y):
    a = 1
    b = 100
    return (a - x)**2 + b * (y - x**2)**2

# データの生成
x = np.linspace(-2, 2, 100)
y = np.linspace(-1, 3, 100)
X, Y = np.meshgrid(x, y)
Z = rosenbrock(X, Y)

# 3Dグラフのプロット
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')

# グラフの表示
plt.show()

このコードでは、ローゼンブロック関数を定義し、その関数を使用して3Dプロットを行っています。

[実行結果]

ソースコード解説

このソースコードは、Pythonを使用してローゼンブロック関数を3Dプロットするものです。

以下は、コードの章立て詳細な説明です：

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

• import numpy as np:
  NumPyライブラリを np としてインポートします。
  NumPyは数値計算を行うための重要なライブラリであり、このコードでは行列操作や数学的関数を使用します。
• import matplotlib.pyplot as plt:
  Matplotlibライブラリから plt としてサブモジュールをインポートします。
  Matplotlibはグラフを描画するためのライブラリであり、このコードではプロットの作成に使用されます。
• from mpl_toolkits.mplot3d import Axes3D:
  Matplotlibの3Dプロット機能を使用するために、Axes3D サブモジュールをインポートします。

2. ローゼンブロック関数の定義：

• def rosenbrock(x, y):
  ローゼンブロック関数を定義します。
  この関数は2つの変数 x と y を取り、それらの値に基づいて関数値を計算します。
  ローゼンブロック関数は、$ ( f(x, y) = (a - x)^2 + b \cdot (y - x^2)^2 ) $の形をしています。
  この関数は、最小値が$ ( x = a ) $および$ ( y = x^2 ) $のときに$ ( 0 ) $になります。

3. データの生成：

• x = np.linspace(-2, 2, 100):
  区間 ([-2, 2]) を100等分した x の値を生成します。
• y = np.linspace(-1, 3, 100):
  区間 ([-1, 3]) を100等分した y の値を生成します。
• X, Y = np.meshgrid(x, y):
  x と y の値からメッシュグリッドを生成します。
  メッシュグリッドは、 x と y の組み合わせに対する関数の値を計算するために使用されます。
• Z = rosenbrock(X, Y):
  ローゼンブロック関数を X と Y のメッシュグリッドに適用し、対応する関数の値を計算します。

4. 3Dグラフのプロット：

• fig = plt.figure():
  新しい図を作成します。
• ax = fig.add_subplot(111, projection='3d'):
  3Dサブプロットを図に追加します。
• ax.plot_surface(X, Y, Z, cmap='viridis'):
  3D曲面をプロットします。
  X、 Y、 Z はそれぞれ x 座標、 y 座標、関数の値であり、 cmap パラメータはカラーマップを指定します。

5. 軸ラベルの設定：

• ax.set_xlabel('X'), ax.set_ylabel('Y'), ax.set_zlabel('Z'):
  X軸、Y軸、Z軸のラベルを設定します。

6. グラフの表示：

• plt.show():
  プロットを表示します。

これにより、ローゼンブロック関数の3Dプロットが生成され、その形状や特性を視覚的に理解することができます。

結果解説

[実行結果]

このグラフは、ローゼンブロック関数$ ( f(x, y) = (a - x)^2 + b \cdot (y - x^2)^2 ) $の3Dプロットです。
この関数は、2つの変数$ ( x ) $と$ ( y ) $の値に依存しています。

• $X軸$と$Y軸$は、それぞれ変数$ ( x ) $と$ ( y ) $の値域を示しています。
  これらの値は、指定した範囲内で等間隔に設定されています。
• $Z軸$は、関数$ ( f(x, y) ) $の値を表しています。
  つまり、Z軸の値は、各$ ( x ) $と$ ( y ) $の組み合わせに対する関数の値です。
• 3Dプロット上の曲面は、関数$ ( f(x, y) ) $の値に対応しています。
  曲面の高さや形状は、$( x ) $と$ ( y ) $の値によって異なります。
  より高い曲面はより大きな関数の値を示し、より低い曲面はより小さな関数の値を示します。
• 曲面の色は、Z値の大きさに基づいて変化します。
  一般的には、Z値が大きいほど明るい色、Z値が小さいほど暗い色になります。

このグラフを通じて、ローゼンブロック関数の形状や極小値を視覚的に理解することができます。

モノレールの運動

モノレールの運動に関する方程式の一例として、モノレール上を移動する列車の位置$ ( x ) $の時間変化を表す単純なモデルを考えます。

ここでは、列車の加速度$ ( a ) $と初速度$ ( v_0 ) $を用いて、次のような運動方程式を考えます。

$$
x(t) = v_0 t + \frac{1}{2} a t^2
$$

Pythonを使用してこの運動方程式を解き、結果をグラフ化します。

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

import numpy as np
import matplotlib.pyplot as plt

# パラメータ
v0 = 10  # 初速度 (m/s)
a = 2    # 加速度 (m/s^2)

# 時間の設定
t = np.linspace(0, 10, 100)  # 0から10までの時間を100分割

# 位置の計算
x = v0 * t + 0.5 * a * t**2

# 結果をプロット
plt.plot(t, x, 'b-', label='Position')
plt.xlabel('Time (s)')
plt.ylabel('Position (m)')
plt.title('Position of Monorail Train Over Time')
plt.legend(loc='best')
plt.grid()
plt.show()

このコードでは、初速度を$ ( v_0 = 10 , \text{m/s} )$、加速度を$ ( a = 2 , \text{m/s}^2 ) $としています。

$0$から$10$までの時間を$100$分割し、その範囲で列車の位置$ ( x ) $を計算しています。

計算結果をプロットしています。

[実行結果]

ソースコード解説

このソースコードは、Pythonを使用してモノレール上を移動する列車の位置を時間の関数として計算し、グラフ化するものです。

以下では、コードの構造と各部分の機能について詳しく説明します。

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

1
2

import numpy as np
import matplotlib.pyplot as plt

• numpy：数値計算用のライブラリ。
  配列操作や数学関数などを提供。
• matplotlib.pyplot：グラフ描画用のライブラリ。
  グラフの作成やカスタマイズが可能。

2. パラメータの設定

1
2

v0 = 10  # 初速度 (m/s)
a = 2    # 加速度 (m/s^2)

• v0：列車の初速度を表すパラメータ。
  単位はメートル毎秒$（m/s）$。
• a：列車の加速度を表すパラメータ。
  単位はメートル毎秒の二乗$（m/s^2）$。

3. 時間の設定

1

t = np.linspace(0, 10, 100)  # 0から10までの時間を100分割

• np.linspace(0, 10, 100)：$0$から$10$までの時間を$100$等分した配列を生成する。
  これにより、$0$秒から$10$秒までの時間を$100$等分して、各時間の値が生成される。

4. 位置の計算

1

x = v0 * t + 0.5 * a * t**2

• x：列車の位置を表す配列。
  初速度と時間の関数として位置を計算し、配列に格納する。
• 位置$ ( x ) $は、初速度$ ( v_0 ) $と時間$ ( t ) $、加速度$ ( a ) $を用いた運動方程式から計算される。

5. 結果のプロット

1
2
3
4
5
6
7

plt.plot(t, x, 'b-', label='Position')
plt.xlabel('Time (s)')
plt.ylabel('Position (m)')
plt.title('Position of Monorail Train Over Time')
plt.legend(loc='best')
plt.grid()
plt.show()

• plt.plot(t, x, 'b-', label='Position')：時間$ ( t ) $に対する位置$ ( x ) $のグラフを作成し、青色の実線でプロットする。ラベルは’Position’として指定される。
• plt.xlabel('Time (s)')：x軸（横軸）のラベルを設定する。
  単位は秒（s）。
• plt.ylabel('Position (m)')：y軸（縦軸）のラベルを設定する。
  単位はメートル（m）。
• plt.title('Position of Monorail Train Over Time')：グラフのタイトルを設定する。
• plt.legend(loc='best')：凡例を表示する。
  ‘best’は自動的に最適な位置に凡例を配置するオプション。
• plt.grid()：グリッド線を表示する。
• plt.show()：グラフを表示する。

これらの手順により、モノレール上を移動する列車の位置を時間の関数として計算し、グラフ化しています。

結果解説

[実行結果]

このグラフは、モノレール上を移動する列車の位置を時間の関数として表しています。
具体的には、縦軸が列車の位置（メートル）、横軸が時間（秒）を表しています。

青い線がグラフ上に描かれており、これは時間とともに列車の位置がどのように変化するかを示しています。
列車は最初の位置から始まり、時間が経つにつれて加速度 $ ( a ) $の影響で移動し、初速度 $ ( v_0 ) $も考慮されます。

グラフが直線の二次関数の形をしていることから、列車の移動が一様に加速されていることがわかります。
初速度$ ( v_0 ) $と加速度$ ( a ) $の値によって、この直線の傾きや曲率が決まります。

超音波

超音波に関する方程式の例として、超音波の進行方向における音圧を表す方程式を考えてみましょう。

超音波の音圧は、物理的な性質や環境によって影響を受けますが、一般的には次のような式で表されます。

$$
P = P_0 \sin(2\pi f t)
$$

ここで、$( P ) $は超音波の音圧、$ ( P_0 ) $は振幅、$ ( f ) $は周波数、$ ( t ) $は時間です。

この方程式を使って、超音波の音圧を計算し、Pythonでグラフ化する方法を示します。

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

import numpy as np
import matplotlib.pyplot as plt

# 定数の設定
P0 = 1.0  # 振幅
f = 1000  # 周波数 (Hz)
t = np.linspace(0, 0.01, 1000)  # 時間 (0から0.01秒まで)

# 超音波の音圧を計算
P = P0 * np.sin(2 * np.pi * f * t)

# グラフの作成
plt.plot(t, P, label='Sound Pressure')
plt.xlabel('Time (s)')
plt.ylabel('Sound Pressure')
plt.title('Sound Pressure of Ultrasonic Wave')
plt.legend()
plt.grid(True)
plt.show()

このコードは、超音波の音圧を計算し、時間に対する音圧のグラフを作成します。

ここでは、振幅を$1.0$、周波数を$1000Hz$、時間を$0$から$0.01$秒まで$1000$点で設定しています。

得られるグラフは、時間の経過とともに音圧が振動する様子を示しています。

[実行結果]

ソースコード解説

このソースコードは、Pythonで超音波の音圧を計算し、グラフ化するためのものです。

以下では、各部分の詳細な説明を行います。

ライブラリのインポート

1
2

import numpy as np
import matplotlib.pyplot as plt

• numpy（略称: np）: 数値計算を行うためのライブラリ。
  このコードでは、数学的な計算や配列操作のために使用されます。
• matplotlib.pyplot（略称: plt）: グラフの描画に使用されるライブラリ。
  このコードでは、グラフの作成やカスタマイズのために使用されます。

定数の設定

1
2
3

P0 = 1.0  # 振幅
f = 1000  # 周波数 (Hz)
t = np.linspace(0, 0.01, 1000)  # 時間 (0から0.01秒まで)

• P0: 超音波の振幅を表す定数。
  このコードでは、$1.0$に設定されています。
• f: 超音波の周波数を表す定数。
  このコードでは、$1000Hz$に設定されています。
• t: 時間の配列を生成するためのNumPyの関数np.linspace()を使用して、$0$から$0.01$秒までの範囲を$1000$個の等間隔に区切った配列を作成しています。
  これは、グラフのX軸（時間軸）に対応します。

超音波の音圧の計算

1

P = P0 * np.sin(2 * np.pi * f * t)

• 超音波の音圧を計算するために、サイン波関数を使用しています。
  周波数fと時間tを用いて、超音波の振幅P0に対する音圧の変化を計算しています。

グラフの作成と表示

1
2
3
4
5
6
7

plt.plot(t, P, label='Sound Pressure')
plt.xlabel('Time (s)')
plt.ylabel('Sound Pressure')
plt.title('Sound Pressure of Ultrasonic Wave')
plt.legend()
plt.grid(True)
plt.show()

• plt.plot(): 時間と超音波の音圧をプロットします。
  tをX軸に、PをY軸に設定しています。
• plt.xlabel(): X軸のラベルを設定します。
  この場合、時間（秒）を表すラベルです。
• plt.ylabel(): Y軸のラベルを設定します。
  この場合、音圧を表すラベルです。
• plt.title(): グラフのタイトルを設定します。
• plt.legend(): グラフに凡例を表示します。
• plt.grid(): グリッド線を表示します。
• plt.show(): グラフを表示します。

これにより、超音波の音圧が時間とともにどのように変化するかを示すグラフが生成されます。

結果解説

[実行結果]

グラフに表示される内容の詳細を説明します。

X軸（横軸）：時間（s）:

X軸は時間を表しており、$0$秒から$0.01$秒までの範囲で$1000$点のデータが表示されています。
この範囲内で、超音波の音圧の振動がどのように変化するかが示されます。

Y軸（縦軸）：音圧:

Y軸は音圧を表しており、グラフ上の各点はその時点での超音波の音圧を示しています。
音圧は単位なしで表されていますが、一般的に振幅を基準にしています。
振幅が大きいほど音圧も大きくなります。

振幅（Amplitude）:

振幅は、音波や超音波の振動の大きさを表します。
グラフ上の波の振れ幅が振幅を示しており、振れ幅が大きいほど音圧も大きくなります。

周波数（Frequency）:

周波数は、音波や超音波の振動が$1$秒間に何回繰り返されるかを表します。
このグラフの例では、周波数が$1000Hz$（$1$秒間に$1000$回の振動）と設定されています。
したがって、グラフ上の振動は$1$秒間に約$1000$回繰り返されます。

グラフ全体から、超音波の音圧が時間の経過とともに振動していることがわかります。
振幅や周波数を調整することで、振動の性質や周期を変えることができます。

超音波に関する方程式の例をあげて、Pythonでといてグラフ化してください。

ケプラーの法則

惑星に関する方程式の例として、惑星の軌道を記述するケプラーの法則を取り上げます。
第一法則によれば、惑星の軌道は楕円軌道であり、その楕円の中心には太陽が位置しています。

ケプラーの第一法則の方程式は楕円の方程式で表されます:

$$
r(\theta) = \frac{a(1 - e^2)}{1 + e \cdot \cos(\theta)}
$$

ここで、$ ( r ) $は太陽からの距離、$ ( \theta ) $は中心からの角度、$ ( a ) $は長半径、$ ( 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

import numpy as np
import matplotlib.pyplot as plt

# ケプラーの法則に基づく関数
def kepler(theta, a, e):
    r = a * (1 - e**2) / (1 + e * np.cos(theta))
    return r

# 角度の範囲を設定
theta = np.linspace(0, 2 * np.pi, 1000)

# 惑星の軌道のパラメータ
a = 1.0  # 長半径
e = 0.5  # 離心率

# 楕円軌道を計算
r = kepler(theta, a, e)

# 楕円をプロット
plt.figure(figsize=(8, 6))
plt.plot(r * np.cos(theta), r * np.sin(theta), label='Orbit')
plt.scatter(0, 0, color='yellow', label='Sun')  # 太陽の位置をプロット
plt.gca().set_aspect('equal', adjustable='box')
plt.title('Kepler\'s First Law')
plt.xlabel('Distance (AU)')
plt.ylabel('Distance (AU)')
plt.legend()
plt.grid(True)
plt.show()

このコードは、指定された長半径と離心率に基づいて楕円軌道を計算し、それをプロットしています。

[実行結果]

ソースコード解説

このソースコードは、Pythonを使用してケプラーの法則に基づいた楕円軌道を計算し、グラフ化するものです。

以下に各部分の詳細を示します。

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

1
2

import numpy as np
import matplotlib.pyplot as plt

• numpyは数値計算を行うためのライブラリで、配列や数学関数などを提供します。
• matplotlib.pyplotはグラフ描画のためのライブラリで、グラフの作成やカスタマイズができます。

2. ケプラーの法則に基づく関数の定義:

1
2
3

def kepler(theta, a, e):
   r = a * (1 - e**2) / (1 + e * np.cos(theta))
   return r

• kepler関数は、与えられた角度$ ( \theta )$、長半径$ ( a )$、離心率$ ( e ) $に基づいて惑星の楕円軌道の距離$ ( r ) $を計算します。

3. 角度の範囲の設定:

1

theta = np.linspace(0, 2 * np.pi, 1000)

• np.linspace関数を使用して、$0$から$ ( 2\pi ) $までの範囲を等間隔で$1000$点に分割した角度の配列を生成します。

4. 惑星の軌道のパラメータの設定:

1
2

a = 1.0  # 長半径
e = 0.5  # 離心率

• 楕円軌道の長半径$ ( a ) $と離心率$ ( e ) $を設定します。

5. 楕円軌道の計算:

1

r = kepler(theta, a, e)

• kepler関数を使用して、与えられた角度、長半径、離心率に基づいて楕円軌道の距離を計算します。

6. 楕円のプロット:

1
2
3
4
5
6
7
8
9
10

plt.figure(figsize=(8, 6))
plt.plot(r * np.cos(theta), r * np.sin(theta), label='Orbit')
plt.scatter(0, 0, color='yellow', label='Sun')  # 太陽の位置をプロット
plt.gca().set_aspect('equal', adjustable='box')
plt.title('Kepler\'s First Law')
plt.xlabel('Distance (AU)')
plt.ylabel('Distance (AU)')
plt.legend()
plt.grid(True)
plt.show()

• plt.figure()で新しい図を作成し、figsizeで図のサイズを設定します。
• plt.plot()で楕円軌道をプロットします。
• plt.scatter()で太陽の位置をプロットします。
• plt.gca().set_aspect('equal', adjustable='box')で軸のアスペクト比を等しく設定します。
• plt.title()、plt.xlabel()、plt.ylabel()でタイトルと軸ラベルを設定します。
• plt.legend()で凡例を表示します。
• plt.grid(True)でグリッドを表示します。
• plt.show()でグラフを表示します。

これにより、ケプラーの法則に基づいた惑星の楕円軌道がプロットされます。

結果解説

[実行結果]

このグラフは、ケプラーの第一法則に基づいて計算された惑星の楕円軌道を示しています。

楕円軌道:

グラフ上の曲線は、惑星が太陽の周りを周回する軌道を表しています。
この軌道は楕円であり、太陽を焦点とする楕円です。
楕円の中心には太陽があります。

太陽:

グラフ上で黄色い点が描かれています。
これは太陽を表しています。
楕円の中心に太陽があることを示しています。

軌道の形状:

楕円の形状は、長半径$ (a) $と離心率$ (e) $の値によって決まります。
長半径が大きいほど楕円はより円に近くなり、離心率が$0$に近いほど楕円は円に近づきます。
離心率が$1$に近いほど楕円はより細長くなります。

軌道の対称性:

楕円は太陽を中心に対称的です。
これは、惑星が太陽の周りを一定の軌道速度で一定の軌道周期で回ることを意味します。

座標軸:

$x軸$と$y軸$はそれぞれ惑星の軌道上の点の$x座標$と$y座標$を表しています。
この座標系は楕円軌道の形状を正しく表現しています。

このグラフは、惑星が楕円軌道を描くことを直感的に理解しやすくするために使用されます。

国際標準大気モデル（International Standard Atmosphere, ISA）

大気圏に関連する方程式の一例として、大気密度の高度依存性をモデル化することが挙げられます。
大気密度は高度によって変化し、これは航空機の飛行や宇宙船の軌道計算などのアプリケーションで重要です。

国際標準大気モデル（International Standard Atmosphere, ISA）は、大気の温度、圧力、密度、速度などの標準的なプロファイルを提供します。
その中で、大気密度の高度依存性は次のような指数関数で表される場合があります:

$$
\rho(h) = \rho_0 \times \exp\left(-\frac{h}{H}\right)
$$

ここで、$[ \rho(h) ] $は高度$ ( h ) $における大気密度、$[ \rho_0 ] $は基準高度$ ( h = 0 ) $における大気密度、$( H ) $はスケールハイト（指数関数の減衰率を決定する定数）です。

これを使って、大気密度の高度依存性を計算し、グラフ化する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

import numpy as np
import matplotlib.pyplot as plt

# 定数
rho0 = 1.225  # 基準高度における大気密度 (kg/m^3)
H = 8000  # スケールハイト (m)

# 高度に対する大気密度の関数
def atmospheric_density(h):
    return rho0 * np.exp(-h / H)

# 高度の範囲
heights = np.linspace(0, 40000, 1000)  # 0から40000メートルの範囲で1000点のデータを生成

# 大気密度の計算
densities = atmospheric_density(heights)

# グラフのプロット
plt.plot(densities, heights)
plt.xlabel('Density (kg/m^3)')
plt.ylabel('Altitude (m)')
plt.title('Atmospheric Density vs Altitude')
plt.grid(True)
plt.gca().invert_yaxis()  # y軸を上向きにする
plt.show()

このコードは、高度に対する大気密度の変化をグラフ化します。

大気密度は高度に応じて指数関数的に減少することがわかります。

[実行結果]

ソースコード解説

このソースコードは、大気圏における高度と大気密度の関係を可視化するためのプログラムです。

以下に、ソースコードの詳細な説明を示します。

インポート

1
2

import numpy as np
import matplotlib.pyplot as plt

• numpy ライブラリは数値計算や配列操作などの機能を提供します。
• matplotlib.pyplot はグラフを描画するためのライブラリです。

定数

1
2

rho0 = 1.225  # 基準高度における大気密度 (kg/m^3)
H = 8000  # スケールハイト (m)

• rho0 は基準高度 (通常、地表) における大気密度を表します。
• H は指数関数の減衰率を決定するスケールハイトを表します。

高度に対する大気密度の関数

1
2

def atmospheric_density(h):
    return rho0 * np.exp(-h / H)

• atmospheric_density 関数は、与えられた高度 h における大気密度を計算します。
• np.exp() は指数関数を計算する numpy ライブラリの関数です。

高度の範囲

1

heights = np.linspace(0, 40000, 1000)  # 0から40000メートルの範囲で1000点のデータを生成

• np.linspace() は、指定された範囲内に等間隔の値を生成します。
  ここでは、$0$から$40000$の間の$1000$個の高度値を生成しています。

大気密度の計算

1

densities = atmospheric_density(heights)

• 先ほど定義した atmospheric_density 関数を使用して、各高度における大気密度を計算します。

グラフのプロット

1
2
3
4
5
6
7

plt.plot(densities, heights)
plt.xlabel('Density (kg/m^3)')
plt.ylabel('Altitude (m)')
plt.title('Atmospheric Density vs Altitude')
plt.grid(True)
plt.gca().invert_yaxis()  # y軸を上向きにする
plt.show()

• plt.plot() は折れ線グラフを描画します。各点の x 座標は大気密度であり、y 座標は高度です。
• plt.xlabel() と plt.ylabel() はそれぞれ x 軸と y 軸のラベルを設定します。
• plt.title() はグラフのタイトルを設定します。
• plt.grid(True) はグリッドを表示します。
• plt.gca().invert_yaxis() は y 軸を反転させ、地表を下方向に表示します。
• plt.show() はグラフを表示します。

結果解説

[実行結果]

このグラフは「高度に対する大気密度の変化」を示しています。

横軸 (x軸):

大気密度 (Density) が表示されています。
単位は$ kg/m^3 $です。
大気密度は地表から高度が上がるにつれて減少します。

縦軸 (y軸):

高度 (Altitude) が表示されています。
単位はメートル$ (m) $です。
地表から上空に向かって高度が増加します。
一般的に、大気密度は高度が上がるにつれて減少しますが、ここでは縦軸が上向きになっているため、地表が上方向、高度が下方向に表示されています。

青線:

大気密度の高度依存性を表しています。
線は、高度が上がるにつれて大気密度が減少することを示しています。
この減少は指数関数的であり、高度が上がるにつれて急速に減少します。

このグラフは、大気中の物体の運動や航空機の飛行、宇宙船の軌道計算など、様々な航空宇宙工学や気象学のアプリケーションにおいて、大気密度がどのように変化するかを理解するのに役立ちます。

リーマンゼータ関数

３次元方程式の一例として、リーマンゼータ関数のモジュラー形式が挙げられます。

リーマンゼータ関数は複素数$ (s) $に対して次のように定義されます。

$$
\zeta(s) = \sum_{n=1}^{\infty} \frac{1}{n^s}
$$

モジュラー形式についての厳密な表現は複雑ですが、その近似式として次のような方程式があります。

$$
f(x, y) = \sum_{n=1}^{\infty} \left( \frac{1}{n^{2 + ix}} + \frac{1}{n^{2 - ix}} \right) \cos\left( \frac{n\pi y}{2} \right)
$$

これを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

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

# リーマンゼータ関数のモジュラー形式の近似式
def zeta_function(x, y, terms=100):
    result = 0
    for n in range(1, terms + 1):
        result += (1 / (n ** (2 + 1j * x)) + 1 / (n ** (2 - 1j * x))) * np.cos(n * np.pi * y / 2)
    return np.real(result)

# x, y の値の範囲
x = np.linspace(-10, 10, 100)
y = np.linspace(-10, 10, 100)
X, Y = np.meshgrid(x, y)

# リーマンゼータ関数のモジュラー形式の値を計算
Z = zeta_function(X, Y)

# 3Dグラフの描画
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('Zeta Function (Modular Form)')
plt.show()

このコードは、リーマンゼータ関数のモジュラー形式の近似式を計算し、3Dグラフで可視化します。

[実行結果]

ソースコード解説

以下にコードの詳細な説明を示します。

1. NumPyとMatplotlibのインポート:

1
2
3

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

• numpy は数値計算を行うためのライブラリです。
• matplotlib.pyplot はグラフの描画に使用されるライブラリです。
• mpl_toolkits.mplot3d は3Dグラフの描画に必要なツールキットです。

2. リーマンゼータ関数のモジュラー形式の近似式:

1
2
3
4
5

def zeta_function(x, y, terms=100):
   result = 0
   for n in range(1, terms + 1):
       result += (1 / (n ** (2 + 1j * x)) + 1 / (n ** (2 - 1j * x))) * np.cos(n * np.pi * y / 2)
   return np.real(result)

• zeta_function 関数は、リーマンゼータ関数のモジュラー形式の近似式を計算します。
• 近似式は、$(x) $と$ (y) $の２つの変数に依存します。
• 近似式は、与えられた$ (x) $と$ (y) $の値に対して計算され、その結果が result に保存されます。

3. x, y の値の範囲の設定:

1
2
3

x = np.linspace(-10, 10, 100)
y = np.linspace(-10, 10, 100)
X, Y = np.meshgrid(x, y)

• np.linspace() 関数は、指定された範囲内で等間隔の数値配列を生成します。
  ここでは、$-10$から$10$までの範囲を$100$分割した値を生成しています。
• np.meshgrid() 関数は、２つの1次元配列から格子状の座標を生成します。
  ここでは、$x$ と$ y $の値の組み合わせを格納する X と Y を生成しています。

4. リーマンゼータ関数のモジュラー形式の値を計算:

1

Z = zeta_function(X, Y)

• zeta_function 関数を用いて、各点$ (x, y) $におけるリーマンゼータ関数のモジュラー形式の近似値を計算し、Z に保存します。

5. 3Dグラフの描画:

1
2
3
4
5
6
7
8

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('Zeta Function (Modular Form)')
plt.show()

• plt.figure() で新しい図を生成し、fig に保存します。
• fig.add_subplot(111, projection='3d') で3Dサブプロットを生成し、ax に保存します。
• ax.plot_surface() で3Dサーフェスプロットを描画します。ここでは、$ (x) $、$ (y) $、$ (z) $の値が X、Y、Z に格納されていると仮定しています。
• ax.set_xlabel()、ax.set_ylabel()、ax.set_zlabel() でそれぞれ軸ラベルを設定します。
• ax.set_title() でグラフのタイトルを設定します。
• plt.show() でグラフを表示します。

結果解説

[実行結果]

このグラフは、リーマンゼータ関数のモジュラー形式の近似式に基づいて描かれています。
リーマンゼータ関数は数学的に興味深い関数であり、数論や解析数学などの分野で重要な役割を果たしています。

具体的には、グラフは３次元空間上に描かれており、横軸は$(x)$、縦軸は$(y)$となっています。
そして、それぞれの点$(x, y)$における関数の値が、高さ方向に表されています。

リーマンゼータ関数のモジュラー形式は、$(x)$と$(y)$の２つの変数に依存する関数であり、その表現は複雑です。
このグラフでは、$(x)$と$(y)$の値が与えられたときの関数の値が色で表現されています。
色の濃淡は関数の値の大きさを示しており、明るい色ほど関数の値が大きく、暗い色ほど関数の値が小さいことを意味します。

リーマンゼータ関数は数学的に興味深い特性を持ち、特にリーマン予想と呼ばれる未解決問題との関連があります。
このグラフは、そのような数学的な概念を視覚的に理解する手助けとなることができます。

連立微分方程式

4次元方程式の例として、次の連立微分方程式を考えます。

$$
\frac{dx}{dt} = -y
$$
$$
\frac{dy}{dt} = x - z
$$
$$
\frac{dz}{dt} = y - w
$$
$$
\frac{dw}{dt} = z
$$

これは、$(x)$、$(y)$、$(z)$、$(w)$の4つの変数を持つ4次元方程式です。

それでは、この方程式を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

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# 4次元微分方程式
def equations(t, y):
    x, y, z, w = y
    dxdt = -y
    dydt = x - z
    dzdt = y - w
    dwdt = z
    return [dxdt, dydt, dzdt, dwdt]

# 初期条件と時間の設定
t_span = (0, 20)
t_eval = np.linspace(*t_span, 1000)
y0 = [1, 0, 0, 0]

# 方程式を解く
sol = solve_ivp(equations, t_span, y0, t_eval=t_eval)

# 結果をプロット
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(sol.y[0], sol.y[1], sol.y[2], label='Trajectory')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.legend()
plt.show()

このコードは、4次元微分方程式を解き、その解を3次元空間上の軌跡としてグラフにプロットします。

[実行結果]

ソースコード解説

以下に、ソースコードの詳細な説明を示します。

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

1
2
3

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

• NumPyは数値計算を行うためのライブラリ、Matplotlibはグラフ描画を行うためのライブラリ、scipy.integrateモジュールのsolve_ivp関数は常微分方程式を解くために使用されます。

2. 4次元微分方程式の定義

1
2
3
4
5
6
7

def equations(t, y):
    x, y, z, w = y
    dxdt = -y
    dydt = x - z
    dzdt = y - w
    dwdt = z
    return [dxdt, dydt, dzdt, dwdt]

• equations関数では、4つの変数$ (x)$,$ (y)$,$ (z)$,$ (w) $に対する4次元微分方程式が定義されています。

3. 初期条件と時間の設定

1
2
3

t_span = (0, 20)
t_eval = np.linspace(*t_span, 1000)
y0 = [1, 0, 0, 0]

• t_spanでは時間の範囲を指定し、t_evalではその範囲内で解を評価する時間の配列を生成します。
• y0は初期条件を表し、この場合は$ (x = 1)$,$ (y = 0)$,$ (z = 0)$,$ (w = 0) $としています。

4. 方程式を解く

1

sol = solve_ivp(equations, t_span, y0, t_eval=t_eval)

• solve_ivp関数を使用して微分方程式を解きます。
  equations関数で定義された微分方程式、初期条件y0、評価する時間t_evalを引数として渡します。

5. 結果をプロット

1
2
3
4
5
6
7
8

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(sol.y[0], sol.y[1], sol.y[2], label='Trajectory')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.legend()
plt.show()

• 解析結果を3次元グラフにプロットします。
  sol.yは解の配列で、sol.y[0]、sol.y[1]、sol.y[2]はそれぞれ$ (x)$、$(y)$、$(z)$の解を表します。

これを用いて3次元軌跡をプロットし、$x軸$、$y軸$、$z軸$にラベルを付けています。

結果解説

[実行結果]

このグラフは、4次元微分方程式の解を表す3次元の軌跡を示しています。

具体的には、$(x)$、$(y)$、$(z)$の3つの変数が時間に対してどのように変化するかを可視化しています。

• $x軸$は変数$(x)$の値を表し、$y軸$は変数$(y)$の値を表します。
• $z軸$は変数$(z)$の値を表しています。
  しかし、実際には4次元方程式を解いているため、$(w)$変数の値が$z軸$に関連付けられます。
  これは、3次元の空間で4次元方程式の解をプロットするための折衷案です。
• 軌跡は、初期条件から始まり、時間とともに変化する解の経路を示しています。
• 解の軌跡がどのように変化するかによって、4次元方程式の動的な振る舞いを視覚的に理解することができます。

このグラフは、4次元の微分方程式を解析する際に、解の性質や振る舞いを視覚化するのに役立ちます。

数理解析の例題として、オイラー法による2階常微分方程式の解法をPythonで解き、グラフ化することを考えます。

この問題は、物理的シミュレーションや工学の問題など、様々な分野で使用されます。

まず、オイラー法による2階常微分方程式の解法について説明します。
オイラー法は、微分方程式を数値的に解くための方法で、特に時間に依存する物理現象を模擬する際によく使用されます。
この方法では、微分方程式を一連の一階の微分方程式に変換し、それぞれを個別に解いていきます

Pythonでこの問題を解くためには、以下の手順を踏むことになります。

1. 時間の上限$T$と時間刻み幅$h$を定義します。
2. 解を保存するための配列$y1[N+1]$,$ y2[N+1]$,$ t[N+1]$を用意します。
3. 微分方程式を満たす関数$f1(t,y1,y2)$, $f2(t,y1,y2)$を定義します。
4. 初期条件$y1[0] = a1$,$ y2[0] = a2$, $t[0] = 0$を設定します。
5. 時間刻み幅$h$で時間tを進めながら、各時点での解を計算し、配列に保存します。
6. 最後に、matplotlibを使用して計算結果をグラフ化します

以下に、これらの手順を実装した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
37
38

import numpy as np
import matplotlib.pyplot as plt

# 定数と初期条件を設定
T = 10
h = 0.01
a1 = 1
a2 = 0

# 配列を用意
y1 = np.zeros(int(T/h)+1)
y2 = np.zeros(int(T/h)+1)
t = np.zeros(int(T/h)+1)

# 微分方程式を満たす関数を定義
def f1(t, y1, y2):
    return y2

def f2(t, y1, y2):
    return -y1

# 初期条件を設定
y1[0] = a1
y2[0] = a2
t[0] = 0

# 時間を進めながら解を計算
for j in range(1, int(T/h)+1):
    t[j] = j*h
    y1[j] = y1[j-1] + h*f1(t[j-1], y1[j-1], y2[j-1])
    y2[j] = y2[j-1] + h*f2(t[j-1], y1[j-1], y2[j-1])

# 結果をグラフ化
plt.figure()
plt.plot(t, y1, label='y1')
plt.plot(t, y2, label='y2')
plt.legend()
plt.show()

このコードは、オイラー法を用いて2階常微分方程式を解き、結果をグラフ化するものです。

微分方程式の形式や初期条件、時間の上限$T$、時間刻み幅$h$は任意に変更することができます。

[実行結果]

ソースコード解説

このソースコードは、2つの連立常微分方程式を数値的に解いてその解をグラフにプロットするプログラムです。

以下にプログラムの詳細な説明を示します：

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

1
2

import numpy as np
import matplotlib.pyplot as plt

• NumPyとMatplotlibをインポートしています。
  NumPyは数値計算のためのライブラリであり、Matplotlibはグラフ描画のためのライブラリです。

2. 定数と初期条件の設定

1
2
3
4

T = 10
h = 0.01
a1 = 1
a2 = 0

• 解析のための定数や初期条件を設定しています。
  ここでは、解析する時間の範囲を表す$T$、時間の刻み幅を表す$h$、初期値$a1$と$a2$を設定しています。

3. 配列の準備

1
2
3

y1 = np.zeros(int(T/h)+1)
y2 = np.zeros(int(T/h)+1)
t = np.zeros(int(T/h)+1)

• 計算結果を格納するための配列$ y1$,$ y2$,$ t $を用意しています。
  これらの配列は、時間の刻み幅$h$に基づいて計算される結果を保持します。

4. 微分方程式を満たす関数の定義

1
2
3
4
5

def f1(t, y1, y2):
    return y2

def f2(t, y1, y2):
    return -y1

• 解きたい連立常微分方程式を満たす関数$ f1$, $f2 $を定義しています。
  この例では、1つ目の方程式が$dy1/dt = y2$、2つ目の方程式が$dy2/dt = -y1$に対応しています。

5. 初期条件の設定

1
2
3

y1[0] = a1
y2[0] = a2
t[0] = 0

• 初期条件を配列に設定しています。

6. 時間を進めながら解を計算

1
2
3
4

for j in range(1, int(T/h)+1):
    t[j] = j*h
    y1[j] = y1[j-1] + h*f1(t[j-1], y1[j-1], y2[j-1])
    y2[j] = y2[j-1] + h*f2(t[j-1], y1[j-1], y2[j-1])

• Eulerの方法を使用して、連立常微分方程式を数値的に解いています。
  時間を進めながら、前の時間ステップの解を使用して次の解を計算します。

7. 結果のグラフ化

1
2
3
4
5

plt.figure()
plt.plot(t, y1, label='y1')
plt.plot(t, y2, label='y2')
plt.legend()
plt.show()

• 得られた解をMatplotlibを使用してグラフにプロットし、$y1$曲線と$y2$曲線を表示しています。
  $y1$曲線は青色で、$y2$曲線はオレンジ色で表されます。

結果解説

[実行結果]

このグラフは、2つの曲線で構成されています。

1. $ y1 $曲線:

この曲線は、時間に対する$ (y1) $の値を表しています。
最初は初期値$ (a1) $から始まり、時間とともに変化します。
曲線の形状は、微分方程式$ (\frac{dy1}{dt} = y2) $に基づいて決定されます。
$ (y2) $の値が変化すると、$ (y1) $の変化速度も変わります。

2. $ y2 $曲線:

この曲線は、時間に対する$ (y2) $の値を表しています。
同様に、初期値$ (a2) $から始まり、時間とともに変化します。
曲線の形状は、微分方程式$ (\frac{dy2}{dt} = -y1) $に基づいて決定されます。
$ (y1) $の値が変化すると、$ (y2) $の変化速度も変わります。

このグラフを通じて、2つの関数$ (y1) $と$ (y2) $が時間の経過とともにどのように振る舞うかを視覚的に理解できます。