常微分方程式

微分方程式の例として、単純な1階の常微分方程式を取り上げましょう。

以下の例では、微分方程式 $ ( \frac{dy}{dt} = -ky ) $を考えます。

この方程式は指数減衰を表します。

これを解いて、その結果をグラフ化します。

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 odeint

# 微分方程式の定義
def model(y, t, k):
    dydt = -k * y
    return dydt

# 初期条件
y0 = 5.0

# 時間の範囲
t = np.linspace(0, 10, 100)

# パラメータ
k = 0.1

# 微分方程式を解く
solution = odeint(model, y0, t, args=(k,))

# 結果のグラフ化
plt.figure(figsize=(8, 6))
plt.plot(t, solution, label='y(t)')
plt.title('指数減衰：$\\frac{dy}{dt} = -ky$')
plt.xlabel('時間 (t)')
plt.ylabel('y(t)')
plt.legend()
plt.grid(True)
plt.show()

このコードでは、SciPyのodeint関数を使用して微分方程式を解き、結果をグラフ化しています。

[実行結果]

指数減衰が確認できるでしょう。

パラメータ$ (k) $を変更することで解の挙動が変わります。

ソースコード解説

以下はソースコードの詳細な説明です。

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

1
2
3

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

• numpy: 数値計算用ライブラリで、数値配列や数学関数などを提供します。
• matplotlib.pyplot: グラフの描画に使用されるライブラリ。
• scipy.integrate.odeint: 常微分方程式を数値的に解くための関数が含まれているSciPyライブラリの一部。

2. 微分方程式の定義:

1
2
3

def model(y, t, k):
   dydt = -k * y
   return dydt

• model関数は、微分方程式$ ( \frac{dy}{dt} = -ky ) $の右辺を定義しています。
  この関数は odeint 関数によって呼び出されます。

3. 初期条件の設定:

1

y0 = 5.0

• 微分方程式の初期条件$ (y(0)) $を設定します。

4. 時間の範囲の設定:

1

t = np.linspace(0, 10, 100)

• 時間の範囲を$0$から$10$までの区間を$100$点で均等に区切って定義します。

5. パラメータの設定:

1

k = 0.1

• 微分方程式のパラメータ$ (k) $を設定します。

6. 微分方程式を解く:

1

solution = odeint(model, y0, t, args=(k,))

• odeint 関数を使用して微分方程式を解きます。
  初期条件 y0、時間の範囲 t、およびパラメータ k を指定しています。

7. 結果のグラフ化:

1
2
3
4
5
6
7
8

plt.figure(figsize=(8, 6))
plt.plot(t, solution, label='y(t)')
plt.title('指数減衰：$\\frac{dy}{dt} = -ky$')
plt.xlabel('時間 (t)')
plt.ylabel('y(t)')
plt.legend()
plt.grid(True)
plt.show()

• matplotlib を使用して、微分方程式の解 solution を時間に対してプロットし、グラフを表示します。
• グラフのタイトル、軸ラベル、凡例、グリッドが設定されています。

このプログラムは、指数減衰を示す微分方程式を数値的に解き、その解をグラフで視覚化する基本的な例です。

結果解説

[実行結果]

以下はグラフに表示される内容の詳細な説明です。

横軸 (X軸):時間$ (t) $

$0$から$10$の範囲を$100$点で等分しています。

縦軸 (Y軸):変数$ (y(t))$

微分方程式の解であり、指数減衰の挙動を示します。
初期条件$ (y(0)) $は$ 5.0 $としています。

グラフの線:

微分方程式を解いた結果の$ (y(t)) $の変化を表しています。
指数減衰の特徴的な形状で、初めは急速に減少し、後に緩やかな減少に移ります。

グラフのタイトル:指数減衰：$ \frac{dy}{dt} = -ky $

微分方程式の性質と解いたものが指数減衰を示すことを示しています。

X軸およびY軸のラベル:

時間（時間$ (t)$）および変数$ (y(t)) $の説明があります。

凡例 (Legend):$y(t)$

グラフ内の線が何を示しているかを示すラベルです。

グリッド:

グラフ内に格子状の線が描かれ、視認性を向上させています。

このグラフは、微分方程式の解を通して、指数減衰が時間とともにどのように進行するかを示しています。

指数減衰は、物理学や生態学などのさまざまな分野で重要な概念であり、この例ではPythonを用いて簡単に表現されています。

最適化問題

最適化問題の例として、特定の制約条件下での関数の最小化を考えてみましょう。

以下は、2変数の関数の最小化問題です。

この例では、scipy.optimize ライブラリを使用して最適化を行います。

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
39
40
41
42
43
44
45
46

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

# 最小化する関数
def objective_function(x):
    return x[0]**2 + x[1]**2  # 例として x^2 + y^2 を最小化

# 制約条件
constraints = ({'type': 'eq', 'fun': lambda x: x[0] + x[1] - 5})  # 制約条件: x + y = 5

# 初期値
initial_guess = [0, 0]

# 最適化の実行
result = minimize(objective_function, initial_guess, constraints=constraints)

# 結果表示
print("最適解:", result.x)
print("最小値:", result.fun)

# グラフの作成
x_vals = np.linspace(-2, 7, 100)
y_vals = np.linspace(-2, 7, 100)
X, Y = np.meshgrid(x_vals, y_vals)
Z = X**2 + Y**2

# 制約条件を描画
plt.plot(x_vals, 5 - x_vals, label='制約条件: x + y = 5', color='red')

# 最小値を描画
plt.scatter(result.x[0], result.x[1], color='green', marker='x', label='最小値')

# 関数の等高線を描画
contour = plt.contour(X, Y, Z, levels=15, cmap='viridis')
plt.clabel(contour, inline=True, fontsize=8)

# グラフの装飾
plt.title('2変数関数の最小化と制約条件')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)

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

このコードでは、関数$ x^2 + y^2 $を最小化する問題を考え、制約条件として$ x + y = 5 $を設定しています。

最適化結果や最小化される関数のグラフが表示されます。

[実行結果]

異なる関数や制約条件を試すことで、様々な最適化問題に対応できます。

ソースコード解説

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

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

1
2
3

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

• numpy: 数値計算のためのライブラリ。
• scipy.optimize.minimize: 最適化問題を解くための関数。
• matplotlib.pyplot: グラフ描画のためのライブラリ。

2. 最小化する関数の定義:

1
2

def objective_function(x):
    return x[0]**2 + x[1]**2  # 例として x^2 + y^2 を最小化

• x は2要素の配列で、x[0] と x[1] が変数となります。
• $ x[0]^2 + x[1]^2 $を最小化する関数を定義しています。

3. 制約条件の設定:

1

constraints = ({'type': 'eq', 'fun': lambda x: x[0] + x[1] - 5})  # 制約条件: x + y = 5

• 制約条件として、$ x + y = 5 $を設定しています。

4. 初期値の設定:

1

initial_guess = [0, 0]

• 最適化の初期値を [0, 0] に設定します。

5. 最適化の実行:

1

result = minimize(objective_function, initial_guess, constraints=constraints)

• minimize 関数を使用して最適化を実行します。
• objective_function を最小化し、制約条件 constraints の下で最適解を見つけます。

6. 結果の表示:

1
2

print("最適解:", result.x)
print("最小値:", result.fun)

• 最適解と最小値を表示します。

7. グラフの作成:

1
2
3
4

x_vals = np.linspace(-2, 7, 100)
y_vals = np.linspace(-2, 7, 100)
X, Y = np.meshgrid(x_vals, y_vals)
Z = X**2 + Y**2

• グラフ描画用の x, y 座標を生成します。
• Z には$ X^2 + Y^2 $の関数値が格納されます。

8. 制約条件の描画:

1

plt.plot(x_vals, 5 - x_vals, label='制約条件: x + y = 5', color='red')

• 制約条件を赤い線で描画します。

9. 最小値の描画:

1

plt.scatter(result.x[0], result.x[1], color='green', marker='x', label='最小値')

• 最小値の点を緑の “x” マークで描画します。

10. 等高線プロットの描画:

1
2

contour = plt.contour(X, Y, Z, levels=15, cmap='viridis')
plt.clabel(contour, inline=True, fontsize=8)

• 関数の等高線プロットを描画します。

11. グラフの装飾:

1
2
3
4
5

plt.title('2変数関数の最小化と制約条件')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)

• グラフのタイトルや軸ラベル、凡例、グリッドを設定します。

12. グラフの表示:

1

plt.show()

• これまでの設定で作成したグラフを表示します。

このコードは、2変数の関数を最小化し、その結果をグラフに可視化しています。

特に、制約条件下での最小値を探索している点がポイントです。

結果解説

[実行結果]

この最適化問題は、関数$ ( f(x, y) = x^2 + y^2 ) $を最小化するというもので、制約条件として$ ( x + y = 5 ) $が課されています。

最適化アルゴリズムは、この関数を最小化する$ ( x ) と ( y ) $の値を見つけることを試みます。

結果は以下の通りです：

• 最適解: [2.5, 2.5]
• 最小値: 12.5

最適解 $ [2.5, 2.5] $は、制約条件 $ ( x + y = 5 ) $を満たす中で$ ( f(x, y) ) $を最小化する点です。

最小値 $ 12.5 $は、この最小化された関数の値です。

グラフには、等高線プロットが表示されています。

等高線は、関数の高さを示し、最小値がグリーンの “x” マークで示されています。

また、制約条件 $ ( x + y = 5 ) $を赤い線で示しています。

最適解がこの赤い線上にあり、かつ関数を最小化する点であることが確認できます。

この例は非常にシンプルですが、最適化問題の基本的な理解を提供しています。

実際の応用では、より複雑な関数や複数の制約条件がありますが、同様のアプローチが使われます。

2変数の定積分

例題として、2変数の定積分を考えてみましょう。

以下の関数を考えます：

$$
f(x, y) = x^2 + y^2
$$

そして、この関数を$ [x \in [-2, 2], y \in [-2, 2]] $の範囲で積分してみます。

Pythonで積分を計算し、結果をグラフ化するためには、SciPyとmatplotlibを使用します。

以下がそのコード例です：

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 import integrate
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 積分する関数
def integrand(x, y):
    return x**2 + y**2

# 積分範囲
x_range = (-2, 2)
y_range = (-2, 2)

# 積分計算
result, _ = integrate.dblquad(integrand, x_range[0], x_range[1], lambda x: y_range[0], lambda x: y_range[1])

print("積分の結果:", result)

# グラフ化
x_vals = np.linspace(x_range[0], x_range[1], 100)
y_vals = np.linspace(y_range[0], y_range[1], 100)

X, Y = np.meshgrid(x_vals, y_vals)
Z = integrand(X, Y)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, cmap='viridis', rstride=5, cstride=5, alpha=0.7, linewidth=0.5, edgecolors='w')

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('2D Integral Example')

plt.show()

このコードでは、dblquadを使用して2変数の積分を計算しています。

そして、plot_surfaceメソッドを使用して関数の3Dプロットを作成しています。

プログラムを実行すると、積分の結果が表示され、関数がどのように分布しているかがグラフで視覚化されます。

[実行結果]

ソースコード解説

このソースコードは、Pythonを使用して2変数の関数を定積分し、その結果を3Dグラフ化するサンプルコードです。

以下に、コードの各部分の説明を示します。

1
2
3
4

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

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

numpyは数値計算を効率的に行うためのライブラリ、scipyは科学技術計算のためのライブラリ、matplotlibはグラフの描画のためのライブラリ、Axes3Dは3Dグラフを描画するためのサブモジュールです。

1
2
3

# 積分する関数
def integrand(x, y):
    return x**2 + y**2

integrandという関数を定義しています。

この関数は、2変数$ [x, y] $に対して$ [x^2 + y^2] $を返します。

これが積分対象の関数です。

1
2
3

# 積分範囲
x_range = (-2, 2)
y_range = (-2, 2)

x_rangeおよびy_rangeで積分の範囲を指定しています。

この例では、$ [x \in [-2, 2]] $および$ [y \in [-2, 2]] $の範囲で積分を行います。

1
2

# 積分計算
result, _ = integrate.dblquad(integrand, x_range[0], x_range[1], lambda x: y_range[0], lambda x: y_range[1])

integrate.dblquad関数を使用して2重積分を計算しています。
dblquadは2変数関数の定積分を計算する関数で、第一引数に積分対象の関数integrand、次に$ [x] $の範囲、そしてその後に$ [y] $の範囲を指定します。
積分結果はresultに格納されます。

1

print("積分の結果:", result)

計算結果をコンソールに出力しています。

1
2
3
4
5
6

# グラフ化
x_vals = np.linspace(x_range[0], x_range[1], 100)
y_vals = np.linspace(y_range[0], y_range[1], 100)

X, Y = np.meshgrid(x_vals, y_vals)
Z = integrand(X, Y)

グラフを描画するために、$ [x] $と$ [y] $の値を生成しています。
linspace関数を使用して、指定された範囲内に等間隔の点を生成しています。
meshgrid関数を使用して、これらの点から座標行列$ [X, Y] $を作成し、それを用いて関数の値$ [Z] $を計算しています。

1
2
3
4
5
6
7
8
9
10

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, cmap='viridis', rstride=5, cstride=5, alpha=0.7, linewidth=0.5, edgecolors='w')

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('2D Integral Example')

plt.show()

最後に、matplotlibを使用して3Dグラフを描画しています。
plot_surfaceメソッドを使用して、$ [X, Y, Z] $の3次元座標で関数をプロットし、その他のオプションでグラフのスタイルを調整しています。

最終的にplt.show()でグラフを表示します。

結果解説

[実行結果]

積分の結果が$ 42.666666666666664 $と表示されました。

これは、関数$ [ f(x, y) = x^2 + y^2 ] $を範囲$ [x \in [-2, 2], y \in [-2, 2]] $で積分した値です。
この積分は、与えられた領域で関数の値を足し合わせることに相当します。

この場合、積分結果が数値で表示されるため、具体的な意味が直感的には理解しにくいかもしれません。
しかし、この数値は範囲内で関数が取る面積や体積に相当します。

積分が正の値であることから、関数$ [ f(x, y) = x^2 + y^2 ] $が指定された範囲で正の値を取り、範囲内の領域の面積が正であることを示しています。

グラフは、積分対象の関数$ [ f(x, y) = x^2 + y^2 ] $の3Dプロットを示しています。
この関数は放物面であり、中心が原点で、$ [ x ] $と$ [ y ] $の値が大きくなるほど関数の値が増加します。
したがって、グラフは中心から外側に向かって放物線が立体的に広がる形状をしています。

この例での積分とグラフは、数学的な計算と視覚的な理解を結びつけるものであり、積分が関数の形状や領域に対する面積や体積の計算にどのように寄与するかを示しています。

人口動態

人口動態の例題として、シンプルな人口成長モデルを考えてみましょう。

ここでは、離散的な時間ステップでの人口の変化を表すシンプルな差分方程式を使用します。

例として、次のような差分方程式を考えます：

$$
[ P_{t+1} = P_t + r \cdot P_t \cdot (1 - \frac{P_t}{K}) ]
$$

ここで、$ (P_t) $は時刻$ (t) $での人口、$ (r) $は増加率、$ (K) $は環境収容力（持続可能な最大人口）です。

以下に、この人口成長モデルを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 matplotlib.pyplot as plt

def population_growth_model(initial_population, growth_rate, carrying_capacity, generations):
    population = [initial_population]

    for t in range(1, generations):
        next_population = population[-1] + growth_rate * population[-1] * (1 - population[-1] / carrying_capacity)
        population.append(next_population)

    return population

# パラメータ設定
initial_population = 100  # 初期人口
growth_rate = 0.1        # 増加率
carrying_capacity = 500   # 環境収容力
generations = 50          # 世代数

# 人口成長モデルの実行
population_over_time = population_growth_model(initial_population, growth_rate, carrying_capacity, generations)

# 結果のグラフ化
plt.plot(range(generations), population_over_time, marker='o', linestyle='-', color='b')
plt.title('人口成長モデル')
plt.xlabel('世代')
plt.ylabel('人口')
plt.grid(True)
plt.show()

このコードでは、初期人口が$100$人で、増加率が$0.1$、環境収容力が$500$人の場合の人口成長モデルをシミュレートしています。

結果を折れ線グラフで表示しています。

[実行結果]

グラフを通じて、初期急激な成長があり、次第に環境収容力に向かって人口が収束していく様子がわかります。

ソースコード解説

このソースコードは、人口成長モデルを定義し、それを用いて人口の時間変化をシミュレートし、最終的には結果をMatplotlibを使用してグラフ化するものです。

以下に各部分の詳細な説明を提供します。

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

1

import matplotlib.pyplot as plt

Matplotlibのpyplotモジュールをインポートしています。

これはグラフの描画に使用されます。

2. 人口成長モデルの定義

1
2
3
4
5
6
7
8

def population_growth_model(initial_population, growth_rate, carrying_capacity, generations):
    population = [initial_population]

    for t in range(1, generations):
        next_population = population[-1] + growth_rate * population[-1] * (1 - population[-1] / carrying_capacity)
        population.append(next_population)

    return population

population_growth_model関数は、差分方程式に基づくシンプルな人口成長モデルを実装しています。
初期人口から始まり、各世代ごとに人口を計算し、結果をリストとして返します。

3. パラメータ設定

1
2
3
4

initial_population = 100  # 初期人口
growth_rate = 0.1        # 増加率
carrying_capacity = 500   # 環境収容力
generations = 50          # 世代数

モデルのパラメータを設定しています。
初期人口、増加率、環境収容力、およびシミュレーションの世代数が含まれます。

4. 人口成長モデルの実行

1

population_over_time = population_growth_model(initial_population, growth_rate, carrying_capacity, generations)

先に定義した人口成長モデルを実行し、時間経過に伴う人口の変化をリストとして取得します。

5. 結果のグラフ化

1
2
3
4
5
6

plt.plot(range(generations), population_over_time, marker='o', linestyle='-', color='b')
plt.title('人口成長モデル')
plt.xlabel('世代')
plt.ylabel('人口')
plt.grid(True)
plt.show()

Matplotlibを使用して人口の時間変化をグラフに描画します。
横軸は世代数、縦軸は人口数です。
グリッドが表示され、グラフのタイトルと軸ラベルが設定されています。
グラフは点線（'o'）で表示され、色は青（'b'）です。

このコード全体では、シンプルな人口成長モデルを定義し、それを用いて特定のパラメータでの人口変化を可視化しています。

簡単な数理解析

方程式の解を求める簡単な例題を取り上げます。

以下の2次方程式を解く問題を考えましょう。

$$
[ax^2 + bx + c = 0]
$$

ここで、$ (a) $,$ (b) $,$ (c) $は任意の実数です。

解の公式を用いて、この方程式の解を求める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

import cmath  # 複素数の計算が必要な場合に使います

def solve_quadratic_equation(a, b, c):
    # 判別式を計算
    discriminant = b**2 - 4*a*c

    # 判別式が非負なら実数解が存在
    if discriminant >= 0:
        root1 = (-b + discriminant**0.5) / (2*a)
        root2 = (-b - discriminant**0.5) / (2*a)
        return root1, root2
    else:
        # 判別式が負の場合、複素数解が存在
        complex_root1 = (-b + cmath.sqrt(discriminant)) / (2*a)
        complex_root2 = (-b - cmath.sqrt(discriminant)) / (2*a)
        return complex_root1, complex_root2

# 例として、a=1, b=-3, c=2の場合を考える
a = 1
b = -3
c = 2

roots = solve_quadratic_equation(a, b, c)
print("方程式の解:", roots)

このコードでは、解の公式を用いて方程式の解を計算しています。

判別式が非負なら実数解が存在し、判別式が負の場合には複素数解が存在します。

上記の例では、$ (a=1) $,$ (b=-3) $,$ (c=2) $の2次方程式の解を求めています。

[実行結果]

方程式の解: (2.0, 1.0)

ソースコード解説

このソースコードは、2次方程式$ (ax^2 + bx + c = 0) $の解を求めるための関数を定義し、その関数を具体的な例で呼び出しています。

以下に、コードの章立てと詳細な説明を示します。

1. 複素数計算のためのライブラリの導入:

1

import cmath  # 複素数の計算が必要な場合に使います

コードの冒頭で cmath ライブラリを導入しています。
これは、複素数を扱うためのライブラリで、実数解だけでなく複素数解も考慮できるようになります。

2. 2次方程式の解を求める関数の定義:

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

def solve_quadratic_equation(a, b, c):
   # 判別式を計算
   discriminant = b**2 - 4*a*c

   # 判別式が非負なら実数解が存在
   if discriminant >= 0:
       root1 = (-b + discriminant**0.5) / (2*a)
       root2 = (-b - discriminant**0.5) / (2*a)
       return root1, root2
   else:
       # 判別式が負の場合、複素数解が存在
       complex_root1 = (-b + cmath.sqrt(discriminant)) / (2*a)
       complex_root2 = (-b - cmath.sqrt(discriminant)) / (2*a)
       return complex_root1, complex_root2

この関数は、2次方程式の係数 (a), (b), (c) を引数として受け取り、解を返すものです。
判別式を計算し、その値に基づいて解の計算方法を分岐させています。

• 判別式が非負の場合（実数解が存在する場合）：
  実数解を計算して root1 および root2 に代入します。
• 判別式が負の場合（複素数解が存在する場合）：
  複素数解を計算して complex_root1 および complex_root2 に代入します。

最終的に、計算された解が返されます。

3. 例として方程式を解く:

1
2
3
4
5
6
7

# 例として、a=1, b=-3, c=2の場合を考える
a = 1
b = -3
c = 2

roots = solve_quadratic_equation(a, b, c)
print("方程式の解:", roots)

具体的な係数$ (a) $,$ (b) $,$ (c) $を与えて、solve_quadratic_equation 関数を呼び出し、得られた解を出力します。

これにより、与えられた2次方程式の解が実数か複素数かに応じて正しく計算され、結果が出力されます。

仮想通貨価格のモデリング (Geometric Brownian Motion)

一般的な仮想通貨価格のモデリングにおいて、Geometric Brownian Motion (GBM) がよく用いられます。

以下はその例です。

$$
dS = \mu S dt + \sigma S dW
$$

• $ (S) $は仮想通貨の価格
• $ (\mu) $はドリフト（平均収益率）
• $ (\sigma) $はボラティリティ（価格変動の度合い）
• $ (W) $はウィーナープロセス（ブラウン運動）

以下はこのモデルを用いた 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

import numpy as np
import matplotlib.pyplot as plt

def generate_price(mu, sigma, S0, days, dt=0.01):
    np.random.seed(42)
    t = np.arange(0, days, dt)
    n = len(t)
    
    S = np.zeros(n)
    S[0] = S0
    
    for i in range(1, n):
        dW = np.random.normal(0, np.sqrt(dt))
        dS = mu * S[i-1] * dt + sigma * S[i-1] * dW
        S[i] = S[i-1] + dS
    
    return t, S

# パラメータの設定
mu = 0.001  # ドリフト
sigma = 0.02  # ボラティリティ
S0 = 100  # 初期価格
days = 365  # シミュレーション期間

# 価格の生成
t, price_simulation = generate_price(mu, sigma, S0, days)

# グラフの描画
plt.figure(figsize=(10, 6))
plt.plot(t, price_simulation, label='仮想通貨価格')
plt.title('仮想通貨価格のGeometric Brownian Motionモデル')
plt.xlabel('時間')
plt.ylabel('価格')
plt.legend()
plt.show()

このコードでは、GBMモデルに基づいて仮想通貨の価格を生成し、Matplotlibを使用して時間と価格のグラフを描画しています。

パラメータ$ (\mu) $、$ (\sigma) $、$ (S0) $を変更することで異なる結果が得られます。

[実行結果]

ソースコード解説

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

1. 必要なライブラリのインポート

1
2

import numpy as np
import matplotlib.pyplot as plt

• numpyは数値計算のためのライブラリで、行列計算や乱数生成などができます。
• matplotlib.pyplotはグラフの描画に使用されます。

2. 価格生成関数の定義

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

def generate_price(mu, sigma, S0, days, dt=0.01):
    np.random.seed(42)
    t = np.arange(0, days, dt)
    n = len(t)
    
    S = np.zeros(n)
    S[0] = S0
    
    for i in range(1, n):
        dW = np.random.normal(0, np.sqrt(dt))
        dS = mu * S[i-1] * dt + sigma * S[i-1] * dW
        S[i] = S[i-1] + dS
    
    return t, S

• generate_price関数はGBMモデルに基づいて仮想通貨の価格を生成する役割を担います。
• muはドリフト（平均収益率）、sigmaはボラティリティ（価格変動の度合い）、S0は初期価格、daysはシミュレーション期間です。
• dtは時間ステップを表し、デフォルトでは$ 0.01 $に設定されています。
• 生成された価格データと対応する時間軸を返します。

3. パラメータの設定

1
2
3
4

mu = 0.001  # ドリフト
sigma = 0.02  # ボラティリティ
S0 = 100  # 初期価格
days = 365  # シミュレーション期間

• 仮想通貨価格のシミュレーションに使用されるパラメータを設定します。

4. 価格の生成

1

t, price_simulation = generate_price(mu, sigma, S0, days)

• 先ほど定義したgenerate_price関数を呼び出し、仮想通貨価格のデータを生成します。

5. グラフの描画

1
2
3
4
5
6
7

plt.figure(figsize=(10, 6))
plt.plot(t, price_simulation, label='仮想通貨価格')
plt.title('仮想通貨価格のGeometric Brownian Motionモデル')
plt.xlabel('時間')
plt.ylabel('価格')
plt.legend()
plt.show()

• matplotlib.pyplotを使用して、生成された仮想通貨価格データをグラフに描画します。
• グラフには時間に対する価格の変動が可視化されます。
• グラフのタイトルや軸ラベル、凡例なども設定されています。

グラフ解説

[実行結果]

このソースコードによって生成されるグラフは、Geometric Brownian Motion（幾何学ブラウン運動）モデルに基づいて仮想通貨の価格変動を模擬したものです。

以下に、グラフの内容を詳しく説明します。

グラフの構成要素:

1. X軸（時間）:

• グラフの横軸は時間を表しています。
• daysパラメータで指定されたシミュレーション期間内の時間軸がプロットされます。

2. Y軸（価格）:

• グラフの縦軸は仮想通貨の価格を表しています。
• 価格の単位は、モデル内での相対的な値です。

3. 曲線（仮想通貨価格）:

• プロットされた曲線は、Geometric Brownian Motionモデルに基づいて生成された仮想通貨の価格変動を示しています。
• ブラウン運動の性質に基づき、確率的な変動が取り入れられており、連続的かつランダムな変動が観察されます。

4. タイトル:

• グラフのタイトルは「仮想通貨価格のGeometric Brownian Motionモデル」です。
• 使用されている価格変動モデルがGBMであることが示唆されています。

5. 軸ラベル:

• X軸とY軸には、「時間」と「価格」という軸ラベルがそれぞれ設定されています。

6. 凡例:

• グラフには「仮想通貨価格」という凡例が設定されており、プロットされた曲線が何を表しているかを示しています。

内容の解釈:

• グラフは時間の経過に伴う仮想通貨価格の変動を視覚化しています。
• ブラウン運動のランダムな性質により、価格は一定のトレンドを持ちつつも、確率的に変動しています。
• 初期価格$（(S0)）$から始まり、時間が進むにつれて価格が変動していることが観察できます。

このグラフは、GBMモデルに基づく単純な価格変動のシミュレーションであり、実際の市場の複雑な動きを模倣しているわけではありません。

ただし、モデルの理解やパラメータの調整を通じて、価格変動の基本的な性質を観察するのに役立ちます。

ランダムウォークモデル

ランダムウォークモデルは、価格や変数がランダムに変動するモデルで、特に金融市場の価格変動などをモデル化するのによく使用されます。

ランダムウォークの基本的な方程式は以下のように表されます。

$$
[ P_t = P_{t-1} + \epsilon_t ]
$$

ここで、$ ( P_t ) $は時刻$ ( t ) $における価格、$ ( \epsilon_t ) $は平均がゼロのランダムな変数（ノイズ）です。

以下は、このランダムウォークモデルを用いて価格の変動をシミュレートし、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

# ランダムウォークに基づく価格変動のモデル
def random_walk(initial_price, num_steps, volatility):
    prices = [initial_price]
    for _ in range(num_steps):
        price = prices[-1] + np.random.normal(0, volatility)
        prices.append(price)
    return prices

# 初期価格
initial_price = 100

# ステップ数
num_steps = 100

# ボラティリティ（価格変動の程度）
volatility = 1.0

# 価格変動のシミュレーション
price_simulation = random_walk(initial_price, num_steps, volatility)

# グラフを描画
plt.plot(range(num_steps + 1), price_simulation)
plt.title('価格変動のランダムウォーク')
plt.xlabel('ステップ')
plt.ylabel('価格')
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
4
5
6

def random_walk(initial_price, num_steps, volatility):
   prices = [initial_price]
   for _ in range(num_steps):
       price = prices[-1] + np.random.normal(0, volatility)
       prices.append(price)
   return prices

• random_walk関数は、初期価格、ステップ数、ボラティリティを受け取り、ランダムウォークモデルに基づいて価格変動をシミュレーションします。
• np.random.normal(0, volatility)は平均が0で標準偏差がvolatilityの正規分布からランダムな変数を生成しています。

3. パラメータの設定:

1
2
3

initial_price = 100
num_steps = 100
volatility = 1.0

• initial_priceは初期価格、num_stepsはステップ数、volatilityは価格変動の程度を表します。

4. 価格変動のシミュレーション:

1

price_simulation = random_walk(initial_price, num_steps, volatility)

• random_walk関数を使用して価格変動をシミュレーションし、その結果をprice_simulationに格納します。

5. グラフの描画:

1
2
3
4
5
6

plt.plot(range(num_steps + 1), price_simulation)
plt.title('価格変動のランダムウォーク')
plt.xlabel('ステップ')
plt.ylabel('価格')
plt.grid(True)
plt.show()

• plt.plotで価格変動のランダムウォークを描画します。
• plt.titleでグラフにタイトルを追加し、plt.xlabelとplt.ylabelで軸のラベルを指定します。
• plt.grid(True)でグリッドを表示し、plt.show()でグラフを表示します。

このコードは、単純なランダムウォークモデルを使用して価格変動をシミュレーションし、それをグラフで視覚化するサンプルです。

グラフ解説

[実行結果]

上記のグラフは、ランダムウォークモデルに基づく価格変動のシミュレーションを表しています。

以下はグラフの詳細な説明です。

横軸 (X軸):

ステップ数を示しています。
各ステップは時間の経過を表しており、時刻$ ( t ) $が$ 0 $から$ 100 $まで変化します。

縦軸 (Y軸):

価格を示しています。初期価格は$ 100 $からスタートしており、各ステップでランダムな変動が加わり、価格が変動します。

グラフの形状:

各ステップでの価格変動がランダムであるため、グラフはランダムに上下に動きます。
ノイズ（ランダムな変数）が価格に加わり、その結果として価格がランダムウォークする様子が観察されます。

ボラティリティの影響:

ボラティリティが$ 1.0 $と設定されているため、価格変動が比較的大きくなっています。
ボラティリティが増加すると、価格の変動が激しくなります。

このグラフは、ランダムウォークモデルを用いて生成された価格変動のシミュレーションであり、実際の金融市場の価格変動などをモデル化するために使用される手法の一つです。

天体の動きは、通常、楕円軌道を描くケプラーの法則によって記述されます。

以下に、楕円軌道上を動く天体の運動を表す方程式と、その方程式をグラフ化するPythonコードの例を示します。

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

• $ ( r(\theta) ) $は楕円軌道上の距離を表します。
• $ ( a ) $は長半径（セミメジャーアクシス）を表します。
• $ ( e ) $は離心率を表します。
• $ ( \theta ) $は楕円軌道上の角度を表します。

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 elliptical_orbit(theta, a, e):
    return (a * (1 - e**2)) / (1 + e * np.cos(theta))

# パラメータの設定
a = 5  # 長半径（セミメジャーアクシス）
e = 0.5  # 離心率

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

# 楕円軌道上の距離を計算
r_values = elliptical_orbit(theta_values, a, e)

# 楕円軌道をグラフ化
plt.polar(theta_values, r_values, label='楕円軌道 (a={}, e={})'.format(a, e))
plt.title('楕円軌道のグラフ')
plt.legend()
plt.show()

このコードでは、ケプラーの法則に基づいた楕円軌道の方程式を用いて、指定された長半径$ (a) $と離心率$ (e) $の楕円軌道を極座標系でグラフ化しています。

[実行結果]

パラメータ$ (a) $と$ (e) $を変更することで、異なる楕円軌道を描画できます。

ソースコード解説

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

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

1
2

import numpy as np
import matplotlib.pyplot as plt

• numpyは数値計算をサポートするライブラリで、数学的な操作に使用されます。
• matplotlib.pyplotはグラフの描画に使用されます。

2. 楕円軌道の方程式の定義:

1
2

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

• elliptical_orbit関数は、楕円軌道の方程式を表しています。
  引数として角度 theta、長半径 a、離心率 e を取ります。

3. パラメータの設定:

1
2

a = 5  # 長半径（セミメジャーアクシス）
e = 0.5  # 離心率

• aは楕円軌道の長半径（セミメジャーアクシス）を指定します。
• eは楕円軌道の離心率を指定します。

4. 角度の範囲を指定:

1

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

• linspace関数は、指定された範囲内で等間隔の数値を生成します。
  ここでは、$0$から$ (2\pi) $までの$1000$個の角度を生成しています。

5. 楕円軌道上の距離を計算:

1

r_values = elliptical_orbit(theta_values, a, e)

• 先ほど定義した elliptical_orbit 関数を使用して、各角度における楕円軌道上の距離を計算します。

6. 楕円軌道をグラフ化:

1
2
3
4

plt.polar(theta_values, r_values, label='楕円軌道 (a={}, e={})'.format(a, e))
plt.title('楕円軌道のグラフ')
plt.legend()
plt.show()

• plt.polar関数は極座標系でのプロットを行います。
  楕円軌道の角度と距離を指定して描画します。
• plt.titleでグラフにタイトルを追加し、plt.legendで凡例を表示します。
• plt.showでグラフを表示します。

このコードは、指定されたパラメータに基づいて楕円軌道を描画し、その形状と特性を可視化します。

グラフ解説

[実行結果]

上記のグラフは、ケプラーの法則に基づく楕円軌道を表しています。

このグラフでは、パラメータとして$ (a = 5) $（長半径が5単位）および$ (e = 0.5) $（離心率が0.5）を使用しています。

グラフの特徴:

1. 楕円の形状:

グラフは楕円の形状を示しており、楕円軌道が中心から離れたり近づいたりすることを表しています。

2. 焦点とセンター:

楕円の中心に焦点があり、この焦点からの距離が離心率によって変化します。
離心率が$0$に近いほど、楕円は円に近くなります。

3. 長半径と短半径:

長半径は楕円の中心から最も遠い点までの距離で、この例では$ (a = 5) $です。
楕円の中心から最も近い点までの距離は$ (a(1 - e)) $になります。

4. 楕円の回転:

グラフは楕円が楕円軌道上を回転していることを示しています。
楕円がどの方向に回転しているかは、楕円軌道上の角度$ ( \theta ) $によって定まります。

これらの要素により、楕円軌道がどのように形成され、どのような特性を持っているかがグラフ上で可視化されています。

関数$ ( f(x) = e^{-x} \cdot \sin(x) ) $の極大値を求め、その結果をグラフ化します。

極大値は、関数の微分がゼロになる点で求められます。

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
39
40
41
42
43
44

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

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

# 関数の微分を定義
def f_prime(x):
    return (-np.exp(-x) * np.sin(x)) + np.exp(-x) * np.cos(x)

# 関数を最小化することで極大値を求める
result = minimize_scalar(lambda x: -f(x), bounds=(0, 4), method='bounded')
max_x = result.x
max_y = f(max_x)

# x の範囲を指定
x_values = np.linspace(0, 4, 100)

# 関数の値を計算
y_values = f(x_values)

# グラフを描画
plt.plot(x_values, y_values, label='f(x) = e^{-x} * sin(x)')
plt.scatter(max_x, max_y, color='red', label='極大値')

# グラフにタイトルとラベルを追加
plt.title('関数のグラフと極大値')
plt.xlabel('x')
plt.ylabel('f(x)')

# 凡例を表示
plt.legend()

# グリッドを表示
plt.grid(True)

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

# 求めた極大値の座標を表示
print("極大値の x 値:", max_x)
print("極大値の y 値:", max_y)

このコードでは、SciPyライブラリのminimize_scalar関数を使用して極大値を求め、Matplotlibを使用して関数のグラフと極大値を描画しています。

[実行結果]

ソースコード解説

以下にソースコードの各部分を詳しく説明します。

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

1
2
3

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

• matplotlib.pyplot: グラフ描画のためのライブラリ。
• numpy: 数値計算のためのライブラリ。
• scipy.optimize.minimize_scalar: 関数の最小化を行うためのSciPyの関数。
  極大値を求める際に利用される。

2. 関数の定義

1
2

def f(x):
    return np.exp(-x) * np.sin(x)

• 関数 f(x) は$ (e^{-x} \cdot \sin(x)) $を返すように定義されています。

3. 関数の微分の定義

1
2

def f_prime(x):
    return (-np.exp(-x) * np.sin(x)) + np.exp(-x) * np.cos(x)

• 関数 f_prime(x) は関数 f(x) の導関数を表しています。

4. 極大値の計算

1
2
3

result = minimize_scalar(lambda x: -f(x), bounds=(0, 4), method='bounded')
max_x = result.x
max_y = f(max_x)

• minimize_scalar 関数を使用して関数 f(x) を最小化し、その結果を result に保存します。
  ここで、lambda x: -f(x) は最大値を求めるために関数を反転させています。
• bounds=(0, 4) は x の範囲を指定しており、極大値を求める際の探索範囲を$0$から$4$に制限しています。
• 最適化の結果から求めた極大値の x 座標を max_x、それに対応する y 座標を max_y に保存します。

5. グラフの描画

1
2
3
4
5
6
7
8
9
10
11

x_values = np.linspace(0, 4, 100)
y_values = f(x_values)

plt.plot(x_values, y_values, label='f(x) = e^{-x} * sin(x)')
plt.scatter(max_x, max_y, color='red', label='極大値')
plt.title('関数のグラフと極大値')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.legend()
plt.grid(True)
plt.show()

• np.linspace(0, 4, 100) は$0$から$4$までの範囲を$100$点で等分した x の値を生成します。
• 関数の値を計算し、その結果を plt.plot でグラフに描画します。
• plt.scatter を使用して極大値を赤い点で示します。
• グラフのタイトルや軸のラベル、凡例、グリッドなどを追加しています。

6. 求めた極大値の座標の表示

1
2

print("極大値の x 値:", max_x)
print("極大値の y 値:", max_y)

• 最後に、求めた極大値の座標を表示しています。

このコードは、数学的な関数の極大値を効果的に求め、その結果をグラフで視覚化するために役立つものです。

結果解説

[実行結果]

求めた極大値の結果は以下の通りです。

• 極大値の x 値: 0.785398102208326
• 極大値の y 値: 0.32239694194483326

これらの値は、関数$ ( f(x) = e^{-x} \cdot \sin(x) ) $の極大値に対応します。

グラフ上では、赤い点が極大値の位置を示しています。

関数の形状と極大値についての説明:

• 関数$ ( f(x) = e^{-x} \cdot \sin(x) ) $は指数関数$ ( e^{-x} ) $と正弦関数$ ( \sin(x) ) $の積です。
• この関数は x 軸に対して周期的であり、指数関数によって急速に減少します。
• 極大値は、正弦関数の振動が指数関数の急激な減少に打ち消される地点で発生します。
  この場合、$ ( x \approx 0.785 ) $で極大値が現れ、その値は$ ( y \approx 0.322 ) $です。
• グラフ上では、関数が x 軸に沿って減少する一方で、局所的な極大値が確認できます。

このようにして、数理解析を用いて関数の極大値を見つけ、それをグラフで視覚化しました。

Weierstrass関数 3Dグラフ化

Weierstrass関数を使用して3Dグラフを描画するサンプルコードです。

Weierstrass関数は、連続だが全ての点で微分可能でない関数であり、複雑な振る舞いを示します。

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

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

# Weierstrass関数の定義
def weierstrass(x, y):
    a = 0.5
    b = 7.0
    k_max = 100
    result = 0.0

    for k in range(k_max):
        result += np.cos(a**k * np.pi * x) * np.cos(b**k * np.pi * y)

    return result

# メッシュグリッドの生成
x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
x, y = np.meshgrid(x, y)

# Weierstrass関数の計算
z = weierstrass(x, y)

# 3Dグラフの描画
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x, y, z, cmap='viridis', edgecolor='k')

# グラフの設定
ax.set_xlabel('X軸')
ax.set_ylabel('Y軸')
ax.set_zlabel('Z軸')
ax.set_title('Weierstrass Function 3D Surface Plot')

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

この例では、Weierstrass関数の形状を視覚化しています。

Weierstrass関数は、周期的でありながらも、微分可能ではないため、興味深い挙動を示します。

[実行結果]

ソースコード解説

ソースコードの各部分の詳細な説明を示します。

1. Weierstrass関数の定義

1
2
3
4
5
6
7
8
9
10

def weierstrass(x, y):
    a = 0.5
    b = 7.0
    k_max = 100
    result = 0.0

    for k in range(k_max):
        result += np.cos(a**k * np.pi * x) * np.cos(b**k * np.pi * y)

    return result

この関数は、Weierstrass関数を計算するためのものです。
関数は2つの引数 x と y を取り、その位置でのWeierstrass関数の値を計算して返します。
Weierstrass関数は、無限級数で表現され、aおよびbといったパラメータによって挙動が制御されます。
この関数は100項までの無限級数を計算しています。

2. メッシュグリッドの生成

1
2
3

x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
x, y = np.meshgrid(x, y)

linspace関数を使用して、-5 から 5 の範囲を等間隔に区切った 100 個の点を生成し、meshgrid関数を使ってメッシュグリッドを作成しています。
これにより、2次元の座標平面ができます。

3. Weierstrass関数の計算

1

z = weierstrass(x, y)

メッシュグリッド上の各点 (x, y) における *:Weierstrass 関数**の値を計算し、z という配列に格納します。

4. 3Dグラフの描画

1
2
3

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

matplotlibライブラリを使用して3Dグラフを描画しています。
plot_surface関数を使用して、メッシュグリッド上の (x, y, z) ポイントでの曲面をプロットします。
カラーマップは’viridis’、エッジカラーは黒に設定されています。

5. グラフの設定

1
2
3
4

ax.set_xlabel('X軸')
ax.set_ylabel('Y軸')
ax.set_zlabel('Z軸')
ax.set_title('Weierstrass Function 3D Surface Plot')

X軸、Y軸、Z軸のラベルおよびグラフのタイトルを設定しています。

6. グラフの表示

1

plt.show()

最後に、作成した3Dグラフを表示します。