最小二乗法（Linear Regression）

Scipyを使用して現実的な問題を解決し、結果をグラフで視覚化する例を示します。

以下は、一般的な問題の一つである最小二乗法（Linear Regression）を用いた簡単なデータ解析の例です。

問題：

あるウェブサイトの広告費用とその広告によって獲得した売上データがあります。

このデータを使用して、広告費用と売上の関係を最小二乗法を使ってモデル化し、将来の広告費用に対する売上を予測したいと思います。

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

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

# サンプルデータ生成
ad_costs = np.array([100, 200, 300, 400, 500])  # 広告費用（千円）
sales = np.array([50, 80, 110, 150, 200])       # 売上（千円）

# 線形モデルを定義（y = mx + c）
def linear_model(x, m, c):
    return m * x + c

# 最小二乗法によるモデルパラメータの推定
params, covariance = curve_fit(linear_model, ad_costs, sales)

# 推定されたモデルパラメータ
m, c = params

# グラフ化
plt.scatter(ad_costs, sales, label='Actual Data')
plt.plot(ad_costs, linear_model(ad_costs, m, c), 'r', label=f'Predicted Model: y = {m:.2f}x + {c:.2f}')
plt.xlabel('Advertising Costs (thousand yen)')
plt.ylabel('Sales (thousand yen)')
plt.legend()
plt.title('Relationship between Advertising Costs and Sales')
plt.grid(True)
plt.show()

# 例えば、広告費用が600の場合の売上を予測
predicted_sales = linear_model(600, m, c)
print(f'広告費用が600の場合の予測売上: {predicted_sales:.2f} 千円')

この例では、Scipyのcurve_fit関数を使用して広告費用と売上の関係を線形モデルでモデル化しました。

そして、グラフを使ってデータとモデルの適合度を可視化し、広告費用が600の場合の売上を予測しました。

これはScipyを使用したデータ解析とグラフ化の基本的な例です。

Scipyはさまざまな科学的な問題を解決するための豊富なツールを提供しており、特定の問題に合わせて適切なツールを選択できます。

流体力学

Scipyを使用して、流体力学の一部として非圧縮のStokes方程式を数値的に解決し、結果をグラフ化する方法を示します。

問題：非圧縮のStokes方程式の解析と結果の可視化

非圧縮のStokes方程式は、流体の静止または低速な流れを記述するために使用されます。

以下は、この方程式を解いて、流体内の速度分布をプロットする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
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import lil_matrix
from scipy.sparse.linalg import spsolve

# グリッドの設定
L = 1.0  # 正方形領域の一辺の長さ
Nx = Ny = 50  # グリッドの分割数
dx = dy = L / (Nx - 1)  # グリッドステップ
x = np.linspace(0, L, Nx)
y = np.linspace(0, L, Ny)
X, Y = np.meshgrid(x, y)

# 境界条件
u_top = np.ones(Nx)
u_bottom = np.zeros(Nx)
u_left = np.zeros(Ny)
u_right = np.zeros(Ny)

# Stokes方程式の数値解法
N = Nx * Ny
A = lil_matrix((N, N))
B = np.zeros(N)

for i in range(1, Nx - 1):
    for j in range(1, Ny - 1):
        row = i * Ny + j
        A[row, row] = -4 / (dx * dy)
        A[row, row - 1] = 1 / (dy ** 2)
        A[row, row + 1] = 1 / (dy ** 2)
        A[row, row - Ny] = 1 / (dx ** 2)
        A[row, row + Ny] = 1 / (dx ** 2)
        B[row] = 0  # 初期化

# 境界条件を適用
A = A.tocsr()
for i in range(Nx):
    row = i * Ny
    A[row, row] = 1
    B[row] = u_top[i]
    A[row + Ny - 1, row + Ny - 1] = 1
    B[row + Ny - 1] = u_bottom[i]

for j in range(Ny):
    row = j
    A[row, row] = 1
    B[row] = u_left[j]
    A[row + (Nx - 1) * Ny, row + (Nx - 1) * Ny] = 1
    B[row + (Nx - 1) * Ny] = u_right[j]

# 方程式を解く
u = spsolve(A, B)
u = u.reshape((Nx, Ny))

# 結果をプロット
plt.figure(figsize=(8, 6))
plt.contourf(X, Y, u, cmap='viridis')
plt.colorbar()
plt.xlabel('X-axis')
plt.ylabel('Y-axis')
plt.title('Velocity Distribution of the Incompressible Stokes Equation')
plt.show()

このコードは、非圧縮のStokes方程式を解いて、正方形領域内の速度分布を可視化します。

Stokes方程式は流体力学の基本的な方程式の1つであり、この例は数値シミュレーションによって流体の挙動を理解するための高度な問題を示しています。

ソースコード解説

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

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

• まず、必要なライブラリをインポートしています。
  • numpy（npとしてエイリアス）: 数値計算ライブラリ
  • matplotlib.pyplot（pltとしてエイリアス）: グラフ描画ライブラリ
  • scipy.sparse.lil_matrix: 疎行列（Sparse Matrix）を表現するためのライブラリの一部
  • scipy.sparse.linalg.spsolve: 疎行列を用いた線形方程式の数値解法を提供するライブラリの一部

2. グリッドの設定:

• 解析領域のサイズ L と、グリッドの分割数 Nx および Ny を設定しています。
• グリッドステップ dx および dy は領域サイズを分割数で割ったものです。
• np.linspace を使用して、x と y の値を生成し、np.meshgrid を使用して2次元座標格子 X と Y を作成します。

3. 境界条件の設定:

• u_top、u_bottom、u_left、および u_right という配列を用意して、境界条件を設定します。
  これらの配列は、領域の上端、下端、左端、右端の速度を表しています。

4. Stokes方程式の数値解法:

• 非圧縮のStokes方程式の数値解法を開始します。
• 疎行列 A とベクトル B を初期化します。
  これらは線形方程式 Ax = B の係数行列と右辺ベクトルです。
• ネストされたループを使用して、各グリッドポイントに対して差分方程式を適用し、A 行列と B ベクトルを構築します。
• このステップでは、有限差分法を使用して方程式を離散化しています。

5. 境界条件の適用:

• A 行列と B ベクトルに境界条件を適用します。
  これにより、領域の境界上の速度制約が反映されます。

6. 方程式の解法:

• A 行列と B ベクトルから線形方程式 Ax = B を解くために spsolve 関数を使用します。
  これにより、速度フィールド u が計算されます。

7. 結果のプロット:

• 最後に、計算された速度フィールド u を可視化します。
• plt.contourf を使用して、速度の分布をカラーマップとしてプロットします。
  速度が色の濃さで表現されています。
• カラーバー、軸ラベル、タイトルなどのプロットの詳細も設定されています。

このコードは、非圧縮のStokes方程式を数値的に解く方法を示し、その解を視覚的に理解するのに役立ちます。
速度場の分布が非圧縮の流体の挙動を示しており、科学や工学の多くの応用分野で使用されています。

グラフ解説

非圧縮のStokes方程式の数値解を示すグラフについて詳しく説明します。

このグラフは、正方形領域内の非圧縮流体の速度分布を表しています。

非圧縮Stokes方程式は、流体の静止または低速な流れを記述するために使用されます。以下はグラフの詳細です：

1. カラーマップ:

• グラフはカラーマップで表示されており、色が速度の大きさを示しています。
  色の濃さが速度の大きさを表しており、濃い色ほど速度が大きいことを示します。

2. X軸とY軸:

• X軸とY軸は正方形領域内の座標を表しており、グラフの広がりを示しています。
  グリッドの分割数がNx=Ny=50であるため、50x50の領域を表示しています。

3. カラーバー:

• グラフの右側にカラーバーが表示されており、色と速度の対応を示しています。
  カラーバーの値が大きいほど速度が高いことを示しています。

4. タイトル:

• グラフのタイトルは「Velocity Distribution of the Incompressible Stokes Equation」となっており、グラフの内容を要約しています。

このグラフは、流体内の速度が非圧縮のStokes方程式に従ってどのように分布するかを示しています。
特に、このグラフは境界条件（境界上の速度制約）に基づいて計算された速度分布を表しています。

カラーマップの色が一様でないことに注目してください。
これは速度が領域内で異なることを示しており、流体の速度勾配を反映しています。

非圧縮Stokes方程式は、静止した流体または低速流れを扱う際に非常に重要であり、このようなシミュレーションは工学、物理学、流体力学などの分野で使用されます。

ロトカ・ヴォルテラ方程式

Scipyを使用して非線形微分方程式を解決し、結果をグラフ化示します。

この問題では、ロトカ・ヴォルテラ方程式を取り上げます。

これは生態学的な相互作用をモデル化するのに使用される非線形微分方程式です。

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

1
2
3

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

次に、ロトカ・ヴォルテラ方程式を定義します。

この方程式は、2つの種の個体数の時間変化を表します。

$$
\frac{dx}{dt} = ax - bxy
$$

$$
\frac{dy}{dt} = -cy + dxy
$$

ここで、$(x)$と$(y)$は2つの異なる種の個体数、$(a)$から$(d)$はパラメータです。

1
2
3
4
5

def lotka_volterra(t, z, a, b, c, d):
    x, y = z
    dxdt = a * x - b * x * y
    dydt = -c * y + d * x * y
    return [dxdt, dydt]

次に、初期条件とシミュレーションの設定を指定します。

1
2
3
4
5
6
7
8
9
10
11

a = 0.1
b = 0.02
c = 0.3
d = 0.01

# 初期値
x0 = 40
y0 = 9

# シミュレーションの時間範囲
t_span = (0, 200)

solve_ivpを使用してロトカ・ヴォルテラ方程式を数値的に解きます。

1

solution = solve_ivp(lotka_volterra, t_span, [x0, y0], args=(a, b, c, d), t_eval=np.linspace(0, 200, 1000))

最後に、2つの種の個体数を時間とともにプロットします。

1
2
3
4
5
6
7
8
9

plt.figure(figsize=(10, 6))
plt.plot(solution.t, solution.y[0], label='Population of Species 1')
plt.plot(solution.t, solution.y[1], label='Population of Species 2')
plt.xlabel('Time')
plt.ylabel('Population')
plt.title('Simulation of Lotka-Volterra Equations')
plt.grid(True)
plt.legend()
plt.show()

このコードを実行すると、2つの種の個体数が時間とともにどのように変化するかがグラフに表示されます。

ロトカ・ヴォルテラ方程式は生態学的な相互作用をモデル化するために使用され、複雑な生態系の挙動を理解するのに役立ちます。

グラフ解析

先ほどのロトカ・ヴォルテラ方程式のシミュレーション結果を示すグラフについて、詳細な説明を提供します。

このグラフは、時間に対する2つの種の個体数を表しています。
ロトカ・ヴォルテラ方程式は、2つの異なる種が相互作用する生態系をモデル化するための数学モデルです。
ここでは、2つの種の個体数が時間とともにどのように変化するかを示しています。

1. 横軸（X軸）:

• 横軸は時間を表しています。シミュレーションの開始から終了までの時間範囲が表示されています。
  この例では0から200までの時間が考慮されています。

2. 縦軸（Y軸）:

• 縦軸は個体数を表しています。
  個体数は2つの異なる種（種1と種2）の個体数を示しています。

3. 線グラフ:

• グラフには2本の線が表示されています。
  • “種1の個体数”（オレンジの線）：この線は種1の個体数の時間変化を示しています。
  • “種2の個体数”（青い線）：この線は種2の個体数の時間変化を示しています。

4. グラフのタイトル:

• グラフのタイトルは「Simulation of Lotka-Volterra Equations(ロトカ・ヴォルテラ方程式のシミュレーション)」となっており、シミュレーションの内容を示しています。

5. 凡例:

• グラフには凡例が表示されており、それぞれの線が何を表しているかを説明しています。
  “種1の個体数”と”種2の個体数”の2つの項目が含まれています。

このグラフから読み取れる情報は、2つの異なる種が相互作用する場合、それらの種の個体数が時間とともにどのように変動するかを示しています。
一方の種の個体数が増加すると、もう一方の種の個体数が減少し、その逆も同様です。
これはロトカ・ヴォルテラ方程式が捕食者と被食者の関係をモデル化するのに適していることを示しています。

グラフの詳細な値や特定の時間点における挙動については、グラフの具体的な数値を確認することで理解できます。
このようなモデルは生態学や生物学の研究において相互作用や種間競争などの重要な概念を理解するのに役立ちます。

ブラウニアン運動

Scipyを使用して、仮想通貨の価格変動をモデル化し、その結果をグラフ化する問題を考えてみましょう。

問題：
仮想通貨の価格変動をブラウニアン運動（ランダムウォーク）モデルでモデル化し、100日間の価格変動をシミュレーションしてグラフ化する。

以下は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

import numpy as np
import matplotlib.pyplot as plt

# ブラウニアン運動のパラメータ
mu = 0.001  # 平均収益率
sigma = 0.02  # 収益率の標準偏差
dt = 1/252  # 1日あたりの期間

# 初期価格
initial_price = 100.0

# シミュレーションの日数
days = 100

# シミュレーションの時間軸を作成
t = np.arange(0, days, dt)

# ブラウニアン運動のシミュレーション
price = np.zeros(len(t))
price[0] = initial_price

for i in range(1, len(t)):
    drift = mu * price[i-1] * dt
    diffusion = sigma * price[i-1] * np.random.normal(0, np.sqrt(dt))
    price[i] = price[i-1] + drift + diffusion

# 価格のグラフをプロット
plt.figure(figsize=(10, 6))
plt.plot(t, price)
plt.title('Brownian Motion Simulation of Cryptocurrency Prices')
plt.xlabel('Days')
plt.ylabel('Price')
plt.grid(True)
plt.show()

このコードは、仮想通貨の価格変動をブラウニアン運動モデルでシミュレーションし、100日間の価格変動をグラフに表示します。

ランダムな価格変動が観察され、その結果がグラフに反映されます。

これは仮想通貨市場の価格変動の一般的な特徴を示すシミュレーションです。

ソースコード解説

このPythonコードは、ブラウニアン運動（ランダムウォーク）モデルを使用して仮想通貨の価格変動をシミュレーションし、結果をグラフ化するためのプログラムです。

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

1. import numpy as npとimport matplotlib.pyplot as plt：

NumPyとMatplotlibをインポートしています。
NumPyは数値計算を行うために使用し、Matplotlibはグラフを描画するために使用します。

2. ブラウニアン運動のパラメータの設定：

• mu：平均収益率。
  この値は価格の平均的な増加または減少を示します。
• sigma：収益率の標準偏差。
  価格変動のばらつき度合いを表します。
• dt：1日あたりの期間。
  この例では1日を取り扱います。

3. initial_price：

初期価格を設定します。
この例では100.0としています。

4. days：

シミュレーションの日数を指定します。この例では100日間の価格変動をシミュレーションします。

5. t：

時間軸を作成します。0からdaysまでの日数を1日あたりの期間dtで区切った配列を作成します。

6. ブラウニアン運動のシミュレーション：

• price：価格の変動を格納するための配列を初期化します。
• ループを使用して、日数ごとに価格を計算します。
  ブラウニアン運動は平均収益率（mu）と収益率の標準偏差（sigma）に基づいて計算され、ランダム性を持ちます。
  driftは平均的な増加または減少を、diffusionはランダムな変動を表します。
  これらを組み合わせて価格を計算し、次の日の価格に反映させます。

7. 価格のグラフのプロット：

• plt.figure()：新しいグラフのフィギュアを作成します。
  サイズは10x6インチです。
• plt.plot(t, price)：時間軸tと価格データpriceを使用して価格のグラフをプロットします。
• plt.title()、plt.xlabel()、plt.ylabel()：グラフにタイトルと軸ラベルを追加します。
• plt.grid(True)：グリッドを表示します。
• plt.show()：グラフを表示します。

このコードは、ブラウニアン運動モデルを使用して仮想通貨の価格変動をシミュレーションし、シミュレーション結果をグラフに可視化するためのものです。

結果のグラフはランダムウォークの特徴を示し、価格の不確実性とランダム性を反映しています。

価格は上下に変動し、初期価格からの変動が観察されます。

グラフ解析

上記のグラフは、ブラウニアン運動モデルを使用してシミュレーションされた仮想通貨の価格変動を示しています。

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

1. 時間軸 (X軸):

グラフのX軸は日数を表しており、シミュレーションの期間を示しています。
この例では100日間の価格変動をシミュレーションしています。

2. 価格 (Y軸):

グラフのY軸は仮想通貨の価格を表しており、シミュレーションの結果に基づいて価格がプロットされています。

3. 初期価格:

シミュレーションの初日での価格は100.0として始まります。
この価格から始まり、ランダムな要因により変動します。

4. ブラウニアン運動:

グラフはランダムウォーク（ブラウニアン運動）を示しており、価格はランダムな上下運動をします。
運動の中には上昇と下降があり、これは市場での価格変動の不確実性を反映しています。

5. 平均収益率と標準偏差:

シミュレーションに使用される平均収益率(mu)と標準偏差(sigma)は、運動の特性を調整します。
平均収益率は価格の平均的な増加または減少を示し、標準偏差は価格変動のばらつき度合いを示します。
これらのパラメータによって、シミュレーションされる価格変動の性質が変わります。

このグラフは、仮想通貨市場の価格変動がランダムであり、予測が難しいことを示しています。

ブラウニアン運動モデルは市場の不確実性やランダム性を捉えるための一般的なモデルであり、実際の仮想通貨市場の価格変動を表現するために使用されることがあります。

価格は上昇と下降を繰り返し、短期的には予測が難しいことが分かります。

振動する弦

振動する弦の問題をScipyを使用して解いてみましょう。

問題の設定:

• 長さLの均一な弦があり、一端を固定してもう一端を自由に振動させます。
• 弦の初期の形状と初速度分布が与えられ、時間とともに弦の振動が発展します。

以下は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
39
40

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

# パラメータ設定
L = 1.0     # 弦の長さ (メートル)
Tension = 100.0  # 弦の張力 (N)
Density = 0.01  # 弦の線密度 (kg/m)
c = np.sqrt(Tension / Density)  # 波速度

# 初期条件
def initial_condition(x):
    if 0.2 <= x <= 0.4:
        return np.sin(10 * np.pi * (x - 0.2))
    else:
        return 0.0

# 弦の振動を数値的に解く
x_span = [0, L]
x_values = np.linspace(0, L, 1000)
initial_state = [initial_condition(x) for x in x_values]

def string_wave(t, state):
    u = state
    dudt = np.zeros_like(u)
    dudt[1:-1] = c**2 * (u[:-2] - 2*u[1:-1] + u[2:])
    return [u[1], *dudt[1:-1], u[-2]]

solution = solve_ivp(string_wave, [0, 2.0], initial_state, t_eval=np.linspace(0, 2.0, 400))

# 結果の可視化
plt.figure(figsize=(10, 6))
for i in range(0, solution.y.shape[1], 40):
    plt.plot(x_values, solution.y[:, i], label=f't = {solution.t[i]:.2f}s')
plt.title('弦の振動')
plt.xlabel('位置 (メートル)')
plt.ylabel('変位')
plt.legend()
plt.grid(True)
plt.show()

このコードは、弦の振動方程式を数値的に解き、弦の振動が時間とともにどのように発展するかを示すグラフを生成します。
初期条件を与えることで、振動パターンが明確に視覚化されます。

Scipyを使用してこのような物理的な問題を解くことができます。

ソースコード解説

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

1. numpyとmatplotlibから必要なモジュールをインポートします。

また、物理的な問題を数値的に解くためにscipy.integrateからsolve_ivpをインポートします。

2. パラメータを設定します。

• Lは弦の長さ（メートル）を示します。
• Tensionは弦の張力（N）を示します。
• Densityは弦の線密度（kg/m）を示します。
• cは波速度を計算するために張力と線密度から計算されます。

3. initial_condition関数は、弦の初期条件を設定します。

特定の範囲（0.2から0.4まで）で正弦波の初期変位を持ち、それ以外の位置では初期変位がゼロです。

4. 弦の振動を数値的に解くための準備を行います。

• x_spanは空間座標xの範囲を示します。
• x_valuesはx座標の値を1000点に離散化します。
• initial_stateはx_valuesの各点における初期変位を計算します。

5. string_wave関数は、時間と状態（弦の変位）を受け取り、時間に対する弦の変位の微分を計算します。

これは弦の振動を記述する偏微分方程式の数値解法です。

6. solve_ivp関数を使用して弦の振動を数値的に解きます。

時間範囲は0から2.0秒までで、400の時間ステップに分割されます。

7. 結果を可視化するためのグラフを作成します。

• ループを使用して時間の異なるスナップショットをプロットし、弦の変位を表示します。
• グラフのタイトル、軸ラベル、凡例、グリッドを設定し、最終的にグラフを表示します。

このコードは、振動する弦の数値的なシミュレーションを行い、振動の進化を時間に沿って可視化するものです。

弦の初期条件や物理的なパラメータを変更することで、さまざまな振動パターンを観察できます。

グラフ解析

生成された弦の振動グラフは、時間に対する変位（弦の形状）を示しています。

以下にグラフの詳細な説明を提供します。

• x軸（位置）:
  グラフの横軸は、弦の長さを0から1メートル（または他の単位）にわたって示しています。
  これは弦の位置を表します。

• y軸（変位）:
  グラフの縦軸は、弦の各位置での変位を示しています。
  変位は振幅として解釈でき、正の値は弦が上に、負の値は弦が下に振動していることを示します。

• カーブ:
  各カーブは異なる時間（t）での弦の形状を表しています。
  t = 0から始まり、時間が経過するにつれて弦の振動が進化しています。
  異なる時間点での弦の形状が示されており、これにより弦が振動し続ける様子が視覚化されています。

• 初期条件:
  初期条件は、弦の形状を決定します。
  この例では、弦の一部が時間0で振動を開始しています。
  具体的には、0.2から0.4までの範囲で弦が正弦波のように振動しています。

このグラフは、時間と空間における物理的な現象、つまり弦の振動を捉えたものです。
時間とともに振動がどのように変化し、初期条件に基づいて振動パターンが形成されるかを可視化しています。

振動の速さ、振幅、および周波数は、弦の特性に基づいて変化します。
グラフは、時間と空間の関係を理解するのに役立ちます。

伝熱問題

伝熱問題を考えます。

具体的には、熱伝導方程式を解いて、熱伝導のプロセスを理解し、結果を温度分布の2Dグラフで表示します。

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
47
48
49
50
51
52

import numpy as np
import matplotlib.pyplot as plt
from scipy import sparse, linalg

# ドメインの設定
length_x = 10.0  # x方向の長さ
length_y = 5.0   # y方向の長さ
nx = 50          # x方向のグリッド数
ny = 25          # y方向のグリッド数
dx = length_x / (nx - 1)
dy = length_y / (ny - 1)

# 初期条件の設定
initial_temperature = 100.0  # 初期温度
boundary_temperature = 0.0   # 境界温度

# タイムステップとシミュレーション時間の設定
nt = 100       # タイムステップ数
dt = 0.01      # タイムステップの間隔
sim_time = nt * dt

# メッシュグリッドの生成
x = np.linspace(0, length_x, nx)
y = np.linspace(0, length_y, ny)
X, Y = np.meshgrid(x, y)

# 初期温度分布の設定
temperature = np.full((ny, nx), initial_temperature)
temperature[0, :] = boundary_temperature
temperature[-1, :] = boundary_temperature
temperature[:, 0] = boundary_temperature
temperature[:, -1] = boundary_temperature

# 伝熱係数
thermal_conductivity = 0.01

# シミュレーションのメインループ
for t in range(nt):
    # 熱伝導方程式の離散化
    dTdx = (temperature[1:-1, 2:] - 2 * temperature[1:-1, 1:-1] + temperature[1:-1, :-2]) / dx**2
    dTdy = (temperature[2:, 1:-1] - 2 * temperature[1:-1, 1:-1] + temperature[:-2, 1:-1]) / dy**2

    temperature[1:-1, 1:-1] += thermal_conductivity * (dTdx + dTdy) * dt

# 結果の表示
plt.figure(figsize=(10, 5))
plt.contourf(X, Y, temperature, cmap='hot')
plt.colorbar(label='Temperature (°C)')
plt.title(f'Temperature Distribution after {sim_time} seconds')
plt.xlabel('X-axis')
plt.ylabel('Y-axis')
plt.show()

このコードは、熱伝導方程式を解き、2D温度分布を計算し、Matplotlibを使用して視覚化します。

この例では、2D領域内での温度分布の時間変化をシミュレートしています。

ソースコード解説

以下に、コードの詳細な説明を提供します：

1. import ステートメント:

• numpy（npとしてエイリアス）: 数値計算を効率的に行うためのPythonライブラリ。
• matplotlib.pyplot（pltとしてエイリアス）: グラフ描画のためのライブラリ。
• scipy.sparse と scipy.linalg: SciPyライブラリの一部で、疎行列や線形代数演算に使用されます。

2. ドメインの設定:

• length_x と length_y: シミュレーション領域のx方向とy方向の長さを定義します。
• nx と ny: x方向とy方向のグリッドの数を定義します。
• dx と dy: x方向とy方向のグリッドの間隔を計算します。

3. 初期条件の設定:

• initial_temperature: シミュレーション領域内の初期温度を設定します。
• boundary_temperature: 境界条件として、領域の境界の温度を設定します。

4. タイムステップとシミュレーション時間の設定:

• nt: タイムステップ数を設定します。
• dt: タイムステップの時間間隔を設定します。
• sim_time: シミュレーションの総時間を計算します。

5. メッシュグリッドの生成:

• numpy.linspace を使用して、x軸とy軸のメッシュグリッドを生成します。
  これにより、2D領域内の位置情報を取得できます。

6. 初期温度分布の設定:

• numpy.full を使用して、初期温度で埋められた2D配列（temperature）を生成します。
  境界条件に従って、領域の境界に対応する部分を設定します。

7. 伝熱係数:

• thermal_conductivity: 熱伝導係数を設定します。
  物質の伝熱性質を示します。

8. シミュレーションのメインループ:

• 時間ステップごとに熱伝導方程式を離散化して解きます。
  これにより、温度分布が時間とともに変化します。
  temperature 配列が更新されます。

9. 結果の表示:

• matplotlib を使用して、シミュレーション結果を可視化します。
  contourf 関数を使用して、2D温度分布をカラーマップで表示します。
  カラーバー、タイトル、軸ラベルなどが追加され、分布の可視化が行われます。

このコードは、熱伝導方程式を数値的に解いて、シミュレーション領域内の温度分布を時間とともに可視化します。
熱伝導のプロセスや温度変化を理解するのに役立ちます。

グラフ解析

下記グラフは、2次元の熱伝導問題のシミュレーション結果を表しています。

グラフの詳細を説明します：

1. カラーマップ:

グラフはカラーマップで表現されており、温度の値に応じて色が異なります。
このカラーマップでは、暖かい温度が明るい色で表され、寒い温度が暗い色で表されています。

2. カラーバー:

グラフの右側にはカラーバーが表示されています。
カラーバーは、色と温度の対応関係を示します。
カラーバーのラベルは「Temperature (°C)」で、グラフ上の各色が示す温度の範囲を示しています。

3. タイトル:

グラフの上部にはタイトルが表示されています。
タイトルにはシミュレーションが実行された時間が記載されており、「Temperature Distribution after [時間] seconds」という形式で表示されます。

4. 軸ラベル:

グラフのX軸とY軸には「X-axis」と「Y-axis」というラベルが表示されています。
これらの軸ラベルは、グラフの座標軸を示し、X軸は水平方向（通常は空間）、Y軸は垂直方向（通常は空間）を表します。

5. 温度分布:

グラフの中央部分には、2D空間内の温度分布が表示されています。
明るい色の領域は高温を示し、暗い色の領域は低温を示します。
この分布は時間経過に伴って変化し、熱伝導プロセスがどのように進行しているかを示しています。

このグラフは、熱伝導方程式に基づいて計算された温度分布を視覚化するものであり、時間経過に伴う温度の変化を把握するのに役立ちます。
暖房や冷却、材料の熱伝導特性の評価など、さまざまな応用分野で使用できる熱伝導のシミュレーションに関連しています。

物体の自由振動

Scipyを使用して、物体の自由振動をシミュレーションし、その結果をグラフ化してみましょう。
自由振動は、質点がばねに固定されており、外力がない状態での振動を表します。

以下は、Scipyを使用して自由振動をシミュレートし、その結果をグラフで表示する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

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

# 自由振動の微分方程式
def spring_mass_damper(y, t, k, m, c):
    x, x_dot = y
    x_dot_dot = -(k/m) * x - (c/m) * x_dot
    return [x_dot, x_dot_dot]

# パラメータ設定
m = 1.0  # 質量 (kg)
k = 10.0  # ばね定数 (N/m)
c = 0.5  # 減衰係数 (N-s/m)

# 初期条件
initial_state = [1.0, 0.0]  # 初期位置 (m) と初期速度 (m/s)

# 時間ステップ
t = np.linspace(0, 10, 1000)  # 0から10秒まで1000ステップ

# 微分方程式を解く
solution = odeint(spring_mass_damper, initial_state, t, args=(k, m, c))

# 位置の変化をプロット
plt.figure(figsize=(10, 4))
plt.plot(t, solution[:, 0])
plt.xlabel('Time (s)')
plt.ylabel('Displacement (m)')
plt.title('Free Vibration of a Spring-Mass-Damper System')
plt.grid(True)
plt.show()

このコードは、ばねに固定された質点の自由振動をシミュレートし、時間に対する位置の変化をグラフに表示します。
振動の減衰を示すグラフが生成されます。

ソースコード解説

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

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

• numpy：数値計算を行うためのライブラリ。
• matplotlib.pyplot：グラフ描画のためのライブラリ。
• scipy.integrate.odeint：常微分方程式を数値的に解くためのSciPyの関数。

2. 自由振動の微分方程式を定義します:

• spring_mass_damper 関数は、引数として現在の状態 y、時間 t、ばね定数 k、質量 m、減衰係数 c を受け取ります。
  この関数は、質点の位置 x と速度 x_dot の二階微分方程式を表しており、質点の加速度 x_dot_dot を計算して返します。これにより、システムの運動を記述します。

3. システムのパラメータを設定します:

• m：質点の質量（単位: kg）。
• k：ばねのばね定数（単位: N/m）。
• c：減衰装置の減衰係数（単位: N-s/m）。

4. 初期条件を設定します:

• initial_state は、初期位置と初期速度を表します。
  初期位置は1.0メートル（m）、初期速度は0.0メートル/秒（m/s）に設定されています。

5. 時間ステップを生成します:

• t は0から10秒までの時間スパンを表し、その間を等間隔で1000のステップに区切ります。
  これにより、シミュレーションの時間範囲と精度が設定されます。

6. 微分方程式を数値的に解きます:

• odeint 関数を使用して微分方程式を解き、質点の位置と速度の時間変化を計算します。
  結果は solution に格納されます。

7. 位置の変化をプロットします:

• plt.plot を使用して、時間 t に対する位置の変化をグラフにプロットします。
• X軸は時間（秒）、Y軸は質点の位置（メートル）を表します。
• グラフのタイトルや軸ラベルを設定し、グリッド線も表示されます。

このコードを実行すると、ばね-質点-減衰装置系の自由振動のシミュレーションが行われ、振動が減衰しながら時間の経過に伴って表示されるグラフが生成されます。

グラフ解析

以下にグラフの詳細な説明を行います。

• 横軸（X軸）は時間（Time）を表しており、0秒から10秒までの時間スパンが表示されています。
• 縦軸（Y軸）は変位（Displacement）を表しており、質点の位置の変化を示します。単位はメートル（m）です。

グラフの特徴:

1. 初期位置:

グラフの最初の点（t=0秒）で、質点の位置は1.0メートルです。
これは、初期条件として設定した値です。

2. 振動:

グラフは振動しており、時間の経過に伴って質点の位置が変化します。
振動の周期性が見られます。

3. 減衰:

振動は減衰しており、振幅が時間とともに減少しています。
これは、減衰係数（c）によるもので、振動が徐々に収束していることを示しています。

4. 振動の周波数:

振動の周期は初期条件やシステムのパラメータに依存しますが、このグラフでは周期的な振動が観察されます。

このグラフは、物理学の基本的な概念である自由振動を示しており、時間と位置の関係を視覚化しています。

時間が経過するにつれて振動が減衰し、最終的に静止状態に収束することが分かります。
これは、ばねと質点のシステムにおける典型的な振動の振る舞いを示しています。

株価データ

株価データの分析と可視化を行います。

以下のコードは、PythonのPandasを使用して特定の株式の過去の株価データを取得し、Matplotlibで結果を折れ線グラフで可視化する例です。

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

import pandas as pd
import matplotlib.pyplot as plt
import yfinance as yf  # Yahoo Financeから株価データを取得するためのライブラリ

# 株価データを取得します。ここではApple Inc. (AAPL) の過去5年間のデータを取得します。
ticker = 'AAPL'  # AAPLはApple Inc.のティッカーシンボル
start_date = '2018-09-01'
end_date = '2023-09-01'

# Yahoo Financeから株価データを取得
data = yf.download(ticker, start=start_date, end=end_date)

# 終値の折れ線グラフを作成します。
plt.figure(figsize=(12, 6))
plt.plot(data['Adj Close'], label=f'{ticker} Closing Price', color='b')
plt.xlabel('Date')
plt.ylabel('Stock Price')
plt.title(f'{ticker} Closing Price Data for the Past 5 Years')
plt.legend()
plt.grid(True)

# グラフを表示します。
plt.show()

このコードは、特定の株式（ここではApple Inc.の株式）の過去5年間の終値データを取得し、折れ線グラフで可視化します。株価の変動や傾向を視覚的に確認できます。

株式のティッカーシンボルを変更することで、他の株式のデータを取得し分析することもできます。

ソースコード解説

以下はコードの詳細な説明です：

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

• pandas (as pd): データの操作と分析を行うためのライブラリ。
• matplotlib.pyplot (as plt): グラフの描画を行うためのライブラリ。
• yfinance (as yf): Yahoo Financeから株価データを取得するためのライブラリ。

2. 株価データを取得するためのパラメータを設定します:

• ticker: 株式のティッカーシンボル（ここでは’APPL’、つまりApple Inc.の株式）を指定します。
• start_date: 取得したいデータの開始日を指定します。
• end_date: 取得したいデータの終了日を指定します。

3. Yahoo Financeから株価データを取得します:

• yf.download() 関数を使用して、指定したティッカーシンボル（’AAPL’）および日付範囲（’start_date’から’end_date’まで）で株価データをダウンロードし、data変数に格納します。

4. 終値の折れ線グラフを作成します:

• plt.figure() 関数を使用して、グラフのサイズを設定します。
• plt.plot() 関数を使用して、data['Adj Close'] から終値データを折れ線グラフとして描画します。labelパラメータでラベルを設定し、colorパラメータで線の色を設定します。
• plt.xlabel() および plt.ylabel() 関数でx軸とy軸のラベルを設定します。
• plt.title() 関数でグラフのタイトルを設定します。
• plt.legend() 関数で凡例（ラベル）を表示します。
• plt.grid(True) でグリッド線を表示します。

5. 最後に、plt.show() を呼び出してグラフを表示します。

このコードを実行すると、指定した期間内のApple Inc.の株価の終値データが折れ線グラフとして表示されます。
このようなグラフを通じて、株価の変動や傾向を視覚的に確認できます。

最小二乗法

scipyを使用して実際の問題である最小二乗法を解決し、結果をグラフ化します。

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

1
2
3

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

次に、データを作成します。ここでは、実際の問題として、ある物体の自由落下の時間と距離の関係を考えます。

1
2

time = np.array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0])
distance = np.array([0.05, 0.20, 0.45, 0.80, 1.25, 1.80, 2.45, 3.20, 4.05, 5.00])

次に、最小二乗法を使用してデータをフィットさせる関数を定義します。

1
2

def linear_func(x, a, b):
    return a * x + b

最小二乗法を実行し、最適なパラメータを取得します。

1
2

params, params_covariance = curve_fit(linear_func, time, distance)
a, b = params

結果をグラフ化します。

1
2
3
4
5
6
7

plt.scatter(time, distance, label='Data')
plt.plot(time, linear_func(time, a, b), label='Linear Fit')
plt.xlabel('Time')
plt.ylabel('Distance')
plt.title('Linear Fit using Least Squares Method')
plt.legend()
plt.show()

これにより、データポイントに最も適合する直線が表示されます。

グラフ解説

このグラフは、時間と距離の関係を示しています。
x軸は時間を表し、y軸は距離を表しています。
データポイントは青い点で表示されており、それぞれの点は実際の観測値を表しています。

また、青い線は最小二乗法によってフィットされた回帰直線を表しています。
この回帰直線は、データポイントに最も適合する直線を表しており、時間と距離の関係を近似的に表現しています。

回帰直線の傾きと切片は、最小二乗法によって求められた最適なパラメータです。
傾きは直線の傾きを表し、切片はy軸との交点を表しています。
これらのパラメータは、時間と距離の関係を数値的に表現するために使用されます。

このグラフを通じて、時間と距離の関係を視覚的に理解することができます。
データポイントが回帰直線に近い位置にある場合、回帰直線がデータによく適合していることを示しています。
逆に、データポイントが回帰直線から離れている場合、回帰直線がデータに適合していないことを示しています。

このようなグラフを使用することで、時間と距離の関係を視覚化し、傾向やパターンを把握することができます。
また、回帰直線を使用することで、未知の時間における距離を予測することも可能です。

データ補間

Scipyを使用してデータの補間（補完）を行い、結果をグラフ化する例をご紹介します。

例として、以下のようなデータセットを考えます。

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
from scipy.interpolate import interp1d

# サンプルデータの作成
x = np.linspace(0, 10, num=11, endpoint=True)
y = np.cos(-x**2/9.0)

# 補間関数の作成
f = interp1d(x, y, kind='cubic')

# 補間結果の計算
x_interp = np.linspace(0, 10, num=100, endpoint=True)
y_interp = f(x_interp)

# グラフのプロット
plt.plot(x, y, 'o', label='Data')
plt.plot(x_interp, y_interp, label='Interpolation')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Data Interpolation')
plt.legend()

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

このコードでは、numpy.linspace()関数を使用して等間隔のデータ点を作成し、それに対応するyの値を計算します。

次に、scipy.interpolate.interp1d()関数を使用してデータの補間関数を作成します。
ここでは、kind='cubic'を指定して3次のスプライン補間を行っています。

最後に、補間関数を使用して新しいxの値に対応するyの値を計算し、元のデータと補間結果をグラフ化して表示します。

この例では、Scipyを使用してデータの補間を行い、補間結果をグラフ化しています。
補間によって、元のデータに存在しないxの値に対応するyの値を推定することができます。

グラフ解説

上記のグラフは、データの補間結果を示しています。

青い点で示されているのは、元のデータセットです。
このデータセットは、0から10までの範囲で等間隔に取られた11個のデータポイントで構成されています。
各データポイントは、関数 $y = cos(-x^2/9.0)$ に基づいて計算されました。

補間関数を使用して、元のデータセットを補間することで、新しいxの値に対応するyの値を推定します。
補間関数は、3次のスプライン補間を使用しています。
補間結果は、補間されたデータポイントを結ぶ滑らかな曲線で表されます。

補間結果は、オレンジ色の線で示されています。
この曲線は、元のデータポイントをなめらかに結んでおり、元のデータセットの間の値を推定しています。
補間によって、元のデータセットに存在しないxの値に対応するyの値を推定することができます。