カルダノの公式（Cardano's formula）

カルダノの公式（Cardano’s formula）は、2次方程式$ ( ax^2 + bx + c = 0 ) $の解の公式です。

この公式は以下のように表されます：

$$
x = \frac{-b \pm \sqrt{b^2 - 4ac}}{2a}
$$

この公式をPythonで使用して、2次方程式の解を計算し、グラフ化する方法を示します。

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

1
2

import numpy as np
import matplotlib.pyplot as plt

次に、2次方程式の解を計算する関数を定義します。

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

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

    if discriminant > 0:
        # 2つの異なる実数解を持つ場合
        x1 = (-b + np.sqrt(discriminant)) / (2*a)
        x2 = (-b - np.sqrt(discriminant)) / (2*a)
        return x1, x2
    elif discriminant == 0:
        # 重解を持つ場合
        x = -b / (2*a)
        return x, x
    else:
        # 実数解を持たない場合
        return None

この関数 solve_quadratic は、2次方程式$ ( ax^2 + bx + c = 0 ) $の係数$ ( a, b, c ) $を受け取り、解を計算して返します。

解が実数の場合は$ ( (x_1, x_2) ) $の形で返し、重解の場合は$ ( (x, x) ) $の形で返します。

実数解が存在しない場合は None を返します。

次に、計算した解をグラフ化します。

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

def plot_quadratic_solution(a, b, c):
    x_values = np.linspace(-10, 10, 400)  # xの範囲を設定
    y_values = a*x_values**2 + b*x_values + c  # 2次方程式のグラフを作成

    plt.figure(figsize=(8, 6))
    plt.plot(x_values, y_values, label=f'{a}x^2 + {b}x + {c} = 0')  # 2次方程式のグラフをプロット

    # 解の計算
    solutions = solve_quadratic(a, b, c)
    if solutions:
        for solution in solutions:
            plt.plot(solution, 0, 'ro')  # 解を赤い点でプロット

    plt.title('Quadratic Equation')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.legend()
    plt.grid(True)
    plt.axhline(0, color='black', linewidth=0.5)
    plt.axvline(0, color='black', linewidth=0.5)
    plt.show()

この関数 plot_quadratic_solution は、2次方程式のグラフを描画します。
2次方程式の解が存在する場合は、解を赤い点でプロットします。

これらの関数を使って、具体的な2次方程式を解いてグラフ化する例を示します。

1
2
3
4
5
6
7

# 例: 2x^2 - 5x + 2 = 0 の解とグラフ化
a = 2
b = -5
c = 2

# 2次方程式の解を計算し、グラフ化
plot_quadratic_solution(a, b, c)

この例では、2次方程式$ ( 2x^2 - 5x + 2 = 0 ) $の解を計算し、その解を含むグラフを描画しています。

他の2次方程式に対しても同様に関数を使用できます。

グラフ解説

[実行結果]

2次方程式$ ( ax^2 + bx + c = 0 ) $の解とそのグラフを詳しく説明します。

1. グラフの描画:

• プロットされた曲線は、2次方程式$ ( y = ax^2 + bx + c ) $のグラフです。
  この曲線は放物線を表し、その形状は係数$ ( a, b, c ) $の値によって決まります。
• 横軸$ ( x ) $は変数$ ( x ) $の値を示し、縦軸$ ( y ) $はそれに対応する2次方程式の値$ ( y ) $を表します。

2. 解の計算:

• 2次方程式$ ( ax^2 + bx + c = 0 ) $の解は、$ ( x ) $の値で曲線が$ ( y = 0 ) $と交わる点です。
  これらの点が解を表します。
• 解は関数 solve_quadratic(a, b, c) を用いて計算されます。
  解が実数解を持つ場合は$ ( (x_1, x_2) ) $の形で表示され、重解の場合は同じ$ ( x ) $の値が二つの解として示されます。

3. 解のプロット:

• 計算された解は赤い点でグラフ上にプロットされます。
  これらの点は、2次方程式の解を示しています。
• 解の点が曲線上にある場合は、それが方程式の解であることを示します。

4. 曲線の特徴:

• 2次方程式の放物線は、係数$ ( a ) $によって開き方が異なります。
  $ ( a > 0 ) $の場合は上に凸の形状（頂点が最小値）、$ ( a < 0 ) $の場合は下に凸の形状（頂点が最大値）を示します。
• 係数$ ( b ) $は曲線の位置を左右にシフトさせ、係数$ ( c ) $は曲線の位置を上下にシフトさせます。

5. グラフの表示:

• グラフは横軸の範囲や解の位置によって、曲線の形状や解の配置が変わります。
• 横軸の範囲は np.linspace() 関数によって設定され、解が存在する範囲や曲線の描画範囲が決まります。

これらの要素を組み合わせることで、2次方程式の解とそのグラフを理解することができます。

グラフは曲線の形状や解の位置を直感的に示し、数学的な関係を視覚的に理解するのに役立ちます。

レナード-ジョーンズポテンシャル

レナード-ジョーンズポテンシャルは、原子間または分子間の相互作用を表現するための一般的なポテンシャルエネルギーモデルです。

このポテンシャルは、2つの粒子（原子や分子）間の距離に依存しており、分子動力学シミュレーションや物理化学で広く使用されています。

レナード-ジョーンズポテンシャル $ ( V(r) ) $は次の式で表されます：

$$
V(r) = 4 \epsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} \right)^6 \right]
$$

ここで、各パラメータの意味は以下の通りです：

• $( r ) $: 分子間の距離（原子核間の最短距離からの距離）
• $( \sigma ) $: レナード-ジョーンズポテンシャルがゼロになる距離であり、分子間の最小距離を表す
• $( \epsilon ) $: ポテンシャルエネルギーの深さを決定する定数（吸引力または斥力の強さを示す）
• $( V(r) ) $: $( r ) $の関数として定義されるポテンシャルエネルギー

このポテンシャルエネルギー関数は、$( r ) $に対して典型的な二重極を持ちます。
距離$ ( r ) $が$ ( \sigma ) $以下の範囲では、ポテンシャルエネルギーは正の値を持ち、原子や分子が斥力を受けます。
一方、距離$ ( r ) $が$ ( \sigma ) $を超えると、ポテンシャルエネルギーは負の値を取り、原子や分子は引力を感じます。

レナード-ジョーンズポテンシャルは、分子間の相互作用を表す一般的な近似式であり、原子や分子の結合エネルギー、結合距離、相互作用の性質を理解するために広く用いられます。

物理学や化学の分野で、分子の構造や性質を研究する際に重要な役割を果たします。

ソースコード

以下は、レナード-ジョーンズポテンシャルを解いてグラフ化する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

import numpy as np
import matplotlib.pyplot as plt

def lennard_jones_potential(r, epsilon, sigma):
    """
    Calculate the Lennard-Jones potential energy for a given distance r.
    
    Parameters:
        r (float or numpy.ndarray): Distance between particles.
        epsilon (float): Depth of the potential energy well.
        sigma (float): Distance at which the potential energy is zero.
        
    Returns:
        float or numpy.ndarray: Potential energy at distance r.
    """
    # Calculate the terms in the LJ potential
    term1 = (sigma / r)**12
    term2 = (sigma / r)**6
    
    # LJ potential energy
    potential_energy = 4 * epsilon * (term1 - term2)
    
    return potential_energy

# Define parameters
epsilon = 1.0  # Depth of the potential energy well
sigma = 1.0    # Distance at which potential energy is zero

# Generate distance values
r_values = np.linspace(0.9, 3.0, 100)  # Range of distances to evaluate

# Calculate LJ potential energies
lj_potentials = lennard_jones_potential(r_values, epsilon, sigma)

# Plot the LJ potential
plt.figure(figsize=(8, 6))
plt.plot(r_values, lj_potentials, label='Lennard-Jones Potential')
plt.xlabel('Distance (r)')
plt.ylabel('Potential Energy (V)')
plt.title('Lennard-Jones Potential Energy')
plt.axhline(0, color='black', linestyle='--', linewidth=0.5)  # Horizontal line at V=0
plt.axvline(sigma, color='red', linestyle='--', label='Sigma (σ)')  # Vertical line at sigma
plt.legend()
plt.grid(True)
plt.show()

このコードでは、lennard_jones_potential 関数で与えられた距離$ ( r ) $に対するレナード-ジョーンズポテンシャルを計算し、numpy.linspace を使用して一定の範囲内で距離$ ( r ) $を生成します。
その後、計算されたポテンシャルエネルギーをグラフ化しています。

グラフでは、距離$ ( r ) $がレナード-ジョーンズポテンシャルの影響を受ける範囲を示すために、ポテンシャルがゼロになる距離$ ( \sigma ) $が赤色の破線で示されています。

[実行結果]

ソースコード解説

ソースコードを詳しく説明します。

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
15
16
17
18
19
20

def lennard_jones_potential(r, epsilon, sigma):
    """
    Calculate the Lennard-Jones potential energy for a given distance r.
    
    Parameters:
        r (float or numpy.ndarray): Distance between particles.
        epsilon (float): Depth of the potential energy well.
        sigma (float): Distance at which the potential energy is zero.
        
    Returns:
        float or numpy.ndarray: Potential energy at distance r.
    """
    # Calculate the terms in the LJ potential
    term1 = (sigma / r)**12
    term2 = (sigma / r)**6
    
    # LJ potential energy
    potential_energy = 4 * epsilon * (term1 - term2)
    
    return potential_energy

• lennard_jones_potential関数は、与えられた距離 r に対するレナード-ジョーンズポテンシャルエネルギーを計算します。
• r は粒子間の距離で、epsilonはポテンシャルエネルギーの深さを表し、sigmaはポテンシャルエネルギーがゼロになる距離を表します。
• term1とterm2はレナード-ジョーンズポテンシャルの計算に使用される項です。
• potential_energyはレナード-ジョーンズポテンシャルの値を計算し、その値を返します。

3. パラメータの定義

1
2

epsilon = 1.0  # Depth of the potential energy well
sigma = 1.0    # Distance at which potential energy is zero

• epsilonとsigmaはレナード-ジョーンズポテンシャルのパラメータであり、それぞれポテンシャルエネルギーの深さとゼロになる距離を表します。

4. 距離の範囲の生成

1

r_values = np.linspace(0.9, 3.0, 100)  # Range of distances to evaluate

• np.linspace関数を使用して、$0.9$から$3.0$までの範囲を等間隔で$100$点に分割した距離の配列 r_values を生成します。

5. レナード-ジョーンズポテンシャルの計算

1

lj_potentials = lennard_jones_potential(r_values, epsilon, sigma)

• 先ほど定義したlennard_jones_potential関数を使用して、各距離 r_values に対するレナード-ジョーンズポテンシャルエネルギーを計算し、lj_potentialsに格納します。

6. プロットの作成と表示

1
2
3
4
5
6
7
8
9
10

plt.figure(figsize=(8, 6))
plt.plot(r_values, lj_potentials, label='Lennard-Jones Potential')
plt.xlabel('Distance (r)')
plt.ylabel('Potential Energy (V)')
plt.title('Lennard-Jones Potential Energy')
plt.axhline(0, color='black', linestyle='--', linewidth=0.5)  # Horizontal line at V=0
plt.axvline(sigma, color='red', linestyle='--', label='Sigma (σ)')  # Vertical line at sigma
plt.legend()
plt.grid(True)
plt.show()

• matplotlib.pyplotを使用してグラフを描画します。
• plt.plotで距離 r_values に対するレナード-ジョーンズポテンシャルエネルギーをプロットします。
• 軸ラベルやタイトルを設定し、水平線と垂直線を追加しています。
• plt.show()でグラフを表示します。

このプログラムは、レナード-ジョーンズポテンシャルを計算し、距離に対するポテンシャルエネルギーの振る舞いを視覚化するためのものです。

グラフ解説

[実行結果]

このグラフは、レナード-ジョーンズポテンシャル（Lennard-Jones potential）を示しています。
レナード-ジョーンズポテンシャルは、原子間または分子間のポテンシャルエネルギーを表す一般的なモデルであり、分子動力学シミュレーションや物理化学で広く使用されています。

以下は、グラフの詳細説明です:

X軸（Distance (r)):

分子間の距離 $( r ) $を表します。
この距離は、原子間の最小距離から離れるにつれて増加します。
単位は任意の長さの単位（例えば、アングストロームなど）です。

Y軸（Potential Energy (V)):

レナード-ジョーンズポテンシャル $( V(r) ) $を表します。
このポテンシャルエネルギーは、分子間の距離 $( r )$ に対して計算され、分子間相互作用のエネルギーを示します。
単位はエネルギーの単位（例えば、エレクトロンボルトなど）です。

Lennard-Jones Potential Curve:

グラフ上の曲線は、レナード-ジョーンズポテンシャル $( V(r) ) $の値を$ ( r ) $の範囲で表します。
曲線は$ ( r ) $が大きくなるにつれてゼロに収束し、$( r ) $が小さいと急激に増加します。

Sigma (σ)の破線:

赤色の破線は、レナード-ジョーンズポテンシャルがゼロになる分子間の最小距離 $( \sigma ) $を示します。
この距離において、ポテンシャルエネルギーは最小値を持ちます。

グラフの特徴:

レナード-ジョーンズポテンシャル曲線は、分子間の距離に応じてポテンシャルエネルギーがどのように変化するかを示します。
$( r ) $が$ ( \sigma ) $以下の範囲では斥力が支配的で$、( r ) $が大きくなると引力が支配的になります。
曲線の形状は、分子間の相互作用に関する基本的な特性を示しています。

このグラフは、物理化学や分子動力学シミュレーションなどで、原子や分子間の相互作用を理解するための重要なツールとして使用されます。

プランクの放射則（Planck's radiation law）

プランクの放射則（Planck’s radiation law）は、物体の放射スペクトルを記述するための物理学の法則です。

この法則は、黒体放射の強度が波長と温度にどのように依存するかを示しています。

以下では、プランクの放射則を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

import numpy as np
import matplotlib.pyplot as plt
from scipy.constants import h, c, k

# プランクの放射則を計算する関数
def planck_radiation(wavelength, temperature):
    numerator = 2 * h * c**2 / wavelength**5
    denominator = np.exp(h * c / (wavelength * k * temperature)) - 1
    return numerator / denominator

# 波長の範囲を定義 (単位: メートル)
wavelengths = np.linspace(1e-9, 3e-6, 1000)  # 1 nm から 3 μm の範囲

# 温度のリストを定義 (単位: ケルビン)
temperatures = [300, 600, 1200]  # 300 K, 600 K, 1200 K の場合を考える

# プランクの放射則を各温度で計算してグラフ化
plt.figure(figsize=(10, 6))
for temp in temperatures:
    intensity = planck_radiation(wavelengths, temp)
    plt.plot(wavelengths * 1e9, intensity, label=f'T = {temp} K')

plt.xlabel('Wavelength (nm)', fontsize=12)
plt.ylabel('Spectral Intensity (W / m^2 / nm / sr)', fontsize=12)
plt.title('Planck\'s Radiation Law', fontsize=14)
plt.legend()
plt.grid(True)
plt.show()

このPythonコードでは、プランクの放射則を計算するための関数 planck_radiation を定義しました。

この関数は、指定された波長と温度に対する放射強度を計算します。

次に、計算した放射強度を異なる温度についてグラフ化します。

wavelengths は波長の範囲を定義し、temperatures は考慮する温度のリストです。
planck_radiation 関数を使用して各温度における放射スペクトルを計算し、結果をグラフにプロットしています。

[実行結果]

グラフは波長に対する放射強度を示しており、異なる温度に対する放射スペクトルの変化を視覚化しています。
この例は、プランクの放射則をPythonで計算してグラフ化する基本的な方法を示しています。

ソースコード解説

以下は、ソースコードの説明をします。

インポートとライブラリの準備

1
2
3

import numpy as np
import matplotlib.pyplot as plt
from scipy.constants import h, c, k

• numpy パッケージは数値計算に使用されます。
• matplotlib.pyplot パッケージはグラフの作成や可視化に使用されます。
• scipy.constants からはプランク定数 (h)、光速度 (c)、ボルツマン定数 (k) を使用します。

プランクの放射則の関数定義

1
2
3
4

def planck_radiation(wavelength, temperature):
    numerator = 2 * h * c**2 / wavelength**5
    denominator = np.exp(h * c / (wavelength * k * temperature)) - 1
    return numerator / denominator

• planck_radiation 関数は、プランクの放射則に基づいて放射強度を計算します。
• wavelength は波長の配列、temperature は温度です。
• numerator は放射則の分子部分を計算し、denominator は分母部分を計算します。
• 関数は計算結果を返します。

波長の範囲の定義

1

wavelengths = np.linspace(1e-9, 3e-6, 1000)

• np.linspace を使用して、$1 nm $から$ 3 μm $の範囲を$ 1000 $等分した波長の配列を作成します。

温度のリストの定義

1

temperatures = [300, 600, 1200]

• グラフ化する各温度のリストを定義します。
  ここでは$ 300 K$、$600 K$、$1200 K $を考慮します。

グラフの作成とプロット

1
2
3
4
5
6
7
8
9
10
11

plt.figure(figsize=(10, 6))
for temp in temperatures:
    intensity = planck_radiation(wavelengths, temp)
    plt.plot(wavelengths * 1e9, intensity, label=f'T = {temp} K')

plt.xlabel('Wavelength (nm)', fontsize=12)
plt.ylabel('Spectral Intensity (W / m^2 / nm / sr)', fontsize=12)
plt.title('Planck\'s Radiation Law', fontsize=14)
plt.legend()
plt.grid(True)
plt.show()

• plt.figure で新しい図を作成し、温度ごとにプランクの放射則を計算してプロットします。
• plot 関数を使用して波長に対する放射強度をプロットします。波長は$ nm $単位で表示されます。
• 軸ラベルやタイトルを設定し、凡例を表示します。
• 最後に show 関数でグラフを表示します。

このプログラムは、プランクの放射則に基づいて異なる温度での放射スペクトルを視覚化し、温度が高いほど短い波長での放射強度が増加することを示しています。

結果解説

[実行結果]

プランクの放射則を計算し、グラフ化することで表示される内容を詳しく説明します。

1. 横軸 (Wavelength - 波長):

• グラフの横軸は波長を表します。
• 単位はナノメートル ($nm$) で、範囲は$1 nm $から$ 3000 nm$ ($3 μm$) までです。
• 波長が短いほど青紫色に近く、長いほど赤色に近い光を示します。

2. 縦軸 (Spectral Intensity - スペクトル強度):

• グラフの縦軸は放射強度を表します。
• 単位はワット毎平方メートル毎ナノメートル毎立体ラジアン ($W/m^2/nm/sr$) です。
• 各波長における単位立体ラジアンあたりの放射強度が示されます。

3. 曲線 (Intensity vs. Wavelength - 強度対波長):

• 各曲線は異なる温度に対する放射スペクトルを示します。
• 温度が高いほど波長による放射強度のピークが高くなり、波長の変化に対する分布が異なります。
• 低温 ($300 K$) の場合はピークが波長の長い側にあり、高温 ($1200 K$) の場合はピークが波長の短い側にあります。

4. タイトルと凡例:

• グラフのタイトルは “Planck’s Radiation Law” で、プランクの放射則を示しています。
• 各曲線に対する凡例は、対応する温度 ($T$) を示しています。

このグラフは、温度が異なる条件下で物体が放射する光の波長ごとの強度分布を示しています。

温度が上昇すると、光のピークの位置が短い波長側に移動し、放射強度が増加することが観察されます。

これにより、物体の温度が光の放射特性にどのように影響を与えるかが可視化されます。

ベルヌーイ方程式

ベルヌーイ方程式は、流体力学における基本的な方程式の一つであり、流れの速さと圧力の関係を表現します。

ベルヌーイ方程式は以下の形で表されます：

$$
P + \frac{1}{2} \rho v^2 + \rho gh = \text{constant}
$$

ここで、各項の意味は次の通りです：

• $(P)$ は流体の圧力
• $(\rho)$ は流体の密度
• $(v)$ は流体の速度
• $(g)$ は重力加速度
• $(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
26
27
28

import numpy as np
import matplotlib.pyplot as plt

# 定数と流体の特性パラメータ
rho = 1000  # 流体の密度 (kg/m^3)
g = 9.81    # 重力加速度 (m/s^2)
h = 10      # 流体の高さ (m)
constant_term = 100000  # 定数項

# 速度の範囲を定義
v_values = np.linspace(0, 10, 100)  # 0から10までの速度 (m/s)

# 圧力を計算する関数
def calculate_pressure(v):
    return constant_term - 0.5 * rho * v**2 - rho * g * h

# 速度と圧力の関係を計算
pressure_values = calculate_pressure(v_values)

# グラフのプロット
plt.figure(figsize=(8, 6))
plt.plot(v_values, pressure_values, label='Pressure vs. Velocity')
plt.xlabel('Velocity (m/s)', fontsize=12)
plt.ylabel('Pressure (Pa)', fontsize=12)
plt.title('Bernoulli Equation: Pressure vs. Velocity', fontsize=14)
plt.grid(True)
plt.legend()
plt.show()

このコードでは、与えられた流体の密度 $(\rho)$、重力加速度 $(g)$、流体の高さ $(h)$、および定数項を使用して、ベルヌーイ方程式に基づいて速度と圧力の関係を計算しています。

そして、計算結果をグラフ化しています。

得られるグラフは、速度と圧力の関係を示しており、流体力学におけるベルヌーイ方程式の応用を視覚的に理解するのに役立ちます。

[実行結果]

ソースコード解説

ソースコードの詳細な説明を行います。

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

1
2

import numpy as np
import matplotlib.pyplot as plt

• numpyは数値計算を効率的に行うためのライブラリであり、npとしてインポートされています。
• matplotlib.pyplotはデータの可視化（グラフ描画）を行うためのライブラリであり、pltとしてインポートされています。

2. 定数と流体の特性パラメータの定義

1
2
3
4

rho = 1000  # 流体の密度 (kg/m^3)
g = 9.81    # 重力加速度 (m/s^2)
h = 10      # 流体の高さ (m)
constant_term = 100000  # 定数項

• rhoは流体の密度 $（kg/m^3）$を表します。
• gは重力加速度 $（m/s^2）$を表します。
• hは流体の高さ $（m）$を表します。
• constant_termはベルヌーイ方程式の定数項であり、圧力を計算する際に使用されます。

3. 速度の範囲を定義

1

v_values = np.linspace(0, 10, 100)  # 0から10までの速度 (m/s)

• np.linspace(0, 10, 100)は$0$から$10$までの速度を等間隔で$100$個のデータポイントに区切って生成します。

4. 圧力を計算する関数の定義

1
2

def calculate_pressure(v):
    return constant_term - 0.5 * rho * v**2 - rho * g * h

• calculate_pressure(v)は、速度 vに対する圧力を計算する関数です。
• ベルヌーイ方程式に基づいて圧力を計算しています。

5. 速度と圧力の関係を計算

1

pressure_values = calculate_pressure(v_values)

• v_valuesに対してcalculate_pressure関数を適用して、速度に対する圧力の値を計算します。

6. グラフのプロット

1
2
3
4
5
6
7
8

plt.figure(figsize=(8, 6))
plt.plot(v_values, pressure_values, label='Pressure vs. Velocity')
plt.xlabel('Velocity (m/s)', fontsize=12)
plt.ylabel('Pressure (Pa)', fontsize=12)
plt.title('Bernoulli Equation: Pressure vs. Velocity', fontsize=14)
plt.grid(True)
plt.legend()
plt.show()

• plt.figure(figsize=(8, 6))で新しい図を作成し、サイズを指定します。
• plt.plot(v_values, pressure_values, label='Pressure vs. Velocity')で速度と圧力の関係をプロットします。
• plt.xlabel('Velocity (m/s)', fontsize=12)で$x$軸（横軸）のラベルを設定します。
• plt.ylabel('Pressure (Pa)', fontsize=12)で$y$軸（縦軸）のラベルを設定します。
• plt.title('Bernoulli Equation: Pressure vs. Velocity', fontsize=14)でグラフのタイトルを設定します。
• plt.grid(True)でグリッド線を表示します。
• plt.legend()で凡例を表示します。
• plt.show()でグラフを表示します。

このプログラムは、ベルヌーイ方程式に基づいて速度と圧力の関係を計算し、matplotlibを使用してグラフ化しています。

計算された圧力値は、速度が増加するにつれてどのように変化するかを視覚化することができます。

結果解説

[実行結果]

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

1. x軸（横軸）：速度 (Velocity)

• グラフの横軸は速度 $（m/s）$を表しています。
• 速度の範囲は$0$から$10$までの値で、$100$個の等間隔なデータポイントがプロットされます。
• この範囲内で流体の速度がどのように変化するかを示します。

2. y軸（縦軸）：圧力 (Pressure)

• グラフの縦軸は圧力（$Pa$：パスカル）を表しています。
• 圧力はベルヌーイ方程式に基づいて計算され、速度の関数としてプロットされます。
• 速度が変化するとともに、流体の圧力がどのように変動するかを示します。

3. グラフのラベルとタイトル

• $x$軸（横軸）のラベルは「Velocity ($m/s$)」として設定されています。
  これは流体の速度を表します。
• $y$軸（縦軸）のラベルは「Pressure ($Pa$)」として設定されています。
  これは流体の圧力を表します。
• グラフのタイトルは「Bernoulli Equation: Pressure vs. Velocity」として設定されています。
  これはベルヌーイ方程式に基づく圧力と速度の関係を示しています。

4. グリッド線と凡例

• グラフにはグリッド線が表示されており、視覚的なガイドとして役立ちます。
• 凡例（Legend）は「Pressure vs. Velocity」として表示され、プロットされたデータの内容を説明しています。

このグラフは、流体の速度が増加すると圧力がどのように変化するかを示しています。

ベルヌーイ方程式は流体力学で重要な法則であり、このグラフを通じて流体の基本的な特性を理解することができます。

円柱

円柱の方程式は一般的に以下のように表されます：

$$
x^2 + y^2 = r^2
$$

ここで、$( r ) $は円柱の半径を表します。

円柱の方程式をPythonでグラフ化するためには、MatplotlibとNumPyを使用します。
具体的なコード例を示します：

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
from mpl_toolkits.mplot3d import Axes3D

# 円柱の半径
r = 3

# メッシュグリッドの作成
theta = np.linspace(0, 2*np.pi, 100)
z = np.linspace(-5, 5, 100)
Theta, Z = np.meshgrid(theta, z)
X = r * np.cos(Theta)
Y = r * np.sin(Theta)

# 3Dプロット
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, color='skyblue', alpha=0.6)

# グラフの設定
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('Cylinder')

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

このコードでは、円柱の方程式 $ ( x^2 + y^2 = r^2 ) $を使用して円柱の表面をプロットしています。

円柱の半径 $( r ) $を変更することで、円柱のサイズを調整することができます。

グラフは3次元で表示され、円柱の構造が視覚化されます。

[実行結果]

ソースコード解説

ソースコードの詳細を説明します。

1. NumPyの利用:

1
2
3
4
5
6
7
8
9

import numpy as np
``
`NumPy`は数値計算やデータ操作に用いられるPythonのライブラリです。
ここでは**円柱の座標**を計算するために使用します。

#### 2. **Matplotlibの利用**:
```python
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

Matplotlibはデータ可視化のためのPythonライブラリであり、mpl_toolkits.mplot3dモジュールは3次元プロットのための機能を提供します。

3. 円柱の半径を設定:

1

r = 3

円柱の半径 $ ( r ) $を設定します。
この値を変更すると、描画される円柱のサイズが変わります。

4. メッシュグリッドの作成:

1
2
3
4
5

theta = np.linspace(0, 2*np.pi, 100)
z = np.linspace(-5, 5, 100)
Theta, Z = np.meshgrid(theta, z)
X = r * np.cos(Theta)
Y = r * np.sin(Theta)

np.linspace関数を使用して角度 $( \theta ) $と高さ $( z ) $の範囲を定義し、np.meshgrid関数でメッシュグリッドを作成します。

そして、円柱の表面の$ ( x ) $座標と$ ( y ) $座標を計算します。

5. 3Dプロットの作成:

1
2
3

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, color='skyblue', alpha=0.6)

plt.figure()で新しい図を作成し、fig.add_subplot(111, projection='3d')で3次元サブプロットを追加します。

ax.plot_surface()で円柱の表面をプロットします。

6. グラフの設定:

1
2
3
4

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

軸ラベルやタイトルなど、グラフの表示設定を行います。

7. グラフの表示:

1

plt.show()

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

これにより、円柱の表面が3次元グラフとして描画されます。
円柱の側面は青色で描かれ、透明度は$0.6$に設定されています。

円柱の底面や側面の形状が視覚的に表現され、円柱の構造がわかりやすくなります。

マリーバインッケ方程式

マリーバインッケ方程式は、3次元空間における代表的な参照集合の1つです。

この方程式は以下のように定義されます。

$$
z = -0.5 * (x^2 + y^2 - 1)^3 - 0.7 * (x^2 + y^2 - 1)^2 - (x^2 - y^2) / 8
$$

この方程式は、$x$、$y$、$z$の3変数で表された非線形の方程式です。

方程式の右辺は3つの項から構成されています。

1. 第1項:$ -0.5 * (x^2 + y^2 - 1)^3$

• この項は、単位球面 $(x^2 + y^2 = 1)$からの距離の3乗に$-0.5$を乗じたものです。
• 球面から離れるほど値が大きくなります。

2. 第2項: $-0.7 * (x^2 + y^2 - 1)^2 $

• この項は、単位球面からの距離の2乗に$-0.7$を乗じたものです。
• 球面から離れるほど値が大きくなりますが、第1項ほど急激ではありません。

3. 第3項: $-(x^2 - y^2)/8$

• この項は、$x^2 - y^2$、つまり$x$軸と$y$軸に沿った楕円の方程式に$-1/8$を乗じたものです。
• これにより、$x=y$の直線に沿って値が小さくなる効果があります。

これら3つの項の組み合わせにより、マリーバインッケ曲面の特徴的な形状が生まれます。

具体的には以下のような特徴があります:

• 中心が最小値$(z=-1.462)$となる深い窪み
• 8つの放射状に伸びる渦巻き状の突起
• $x$軸と$y$軸に関する対称性
• 複雑な3次元形状

このマリーバインッケ曲面は、フラクタル幾何学や力学系の研究において重要な役割を果たしています。

また、この方程式は最適化問題の目的関数としても使われることがあります。

高度な数学的性質を持つユニークな曲面として知られています。

ソースコード

次のコードは、NumPy と Matplotlib を使用してマリーバインッケ方程式の3D曲面をプロットします。

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

import numpy as np
import matplotlib.pyplot as plt

# マリーバインッケ方程式
def mariebyinke(x, y):
    return -0.5 * (x**2 + y**2 - 1)**3 - 0.7 * (x**2 + y**2 - 1)**2 - (x**2 - y**2) / 8

# グリッドデータの生成
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
Z = mariebyinke(X, Y)

# 3D プロット
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='viridis')
ax.set_title('Mariebyinke Surface')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')

plt.show()

1. mariebyinke 関数は、指定された$ x $と$ y $の値から方程式の値を計算します。
2. np.linspace を使って$ x $と$ y $の値の範囲を設定し、np.meshgrid でグリッドデータを生成します。
3. mariebyinke 関数を使って$ Z $値 (方程式の値) を計算します。
4. plt.figure と add_subplot を使って新しい図とサブプロットを作成します。
5. plot_surface を使って3D曲面をプロットします。
   rstride と cstride は表示の密度を制御します。
6. 軸ラベルとタイトルを設定します。
7. plt.show() でプロットを表示します。

このコードを実行すると、マリーバインッケ曲面の3Dプロットが表示されます。

曲面のカラーマップは viridis を使用していますが、必要に応じて変更できます。

[実行結果]

グラフ解説

生成されるグラフには、以下の特徴があります:

1. 形状

• グラフは、複雑な3次元曲面の形状を持っています。
• 曲面は、中央部分が窪んでいる形状をしており、外側に向かって持ち上がっています。
• 曲面には、中心から放射状に伸びる8つの渦巻き状の突起があります。

2. 対称性

• 曲面は、$x$軸と$y$軸に関して対称な形状をしています。
  つまり、$x=0$と$y=0$を回転軸として、曲面を$180$度回転させると同じ形状になります。

3. 極値

• 曲面の中心$(x=0, y=0)$が最小値となっています。
• 8つの渦巻き状の突起の先端が局所的な極大値となっています。

4. カラーマップ

• デフォルトでは、曲面の高さ($z$値)に応じて、viridisカラーマップが適用されています。
• 青系の色が低い値、赤系の色が高い値を表しています。

5. 軸

• $x$軸、$y$軸、$z$軸がそれぞれ表示されています。
• $x$軸と$y$軸の範囲は$-2$から$2$までの値が表示されています。

6. タイトル

• グラフのタイトルは「Mariebyinke Surface」と表示されています。

7. 視点

• 3D空間上の視点は調整可能で、様々な角度から曲面を観察できます。

このグラフは、マリーバインッケ方程式によって生成された特徴的な3次元曲面の形状を可視化しています。

この曲面は数学的に興味深い性質を持っており、様々な分野で研究されています。

楕円錐の方程式

楕円錐の方程式は次のように表されます:

$$
(z - w)^2 / c^2 = (x - u)^2 / a^2 + (y - v)^2 / b^2
$$

• $(x, y, z) $は空間内の任意の点の座標
• $(u, v, w) $は楕円錐の頂点の座標
• $a$, $b$, $c$ はそれぞれ$x$, $y$, $z$軸方向の半径

この方程式は以下のようにも書くことができます:

$$
z = w ± c/a * sqrt(a^2 - (x - u)^2 - (y - v)^2)
$$

符号の$+$ は$z >= w$の上半分、符号の$-$ は$z <= w$の下半分を表します。

この方程式は、$(x, y)$平面上の点$(x-u, y-v)$から頂点$(u, v, w)$を通る直線の長さをcで割ったものが、楕円錐の半径 $a$, $b$の比になることを表しています。

つまり、楕円錐とは、頂点から$x-y$平面に下ろした垂線の長さと、その垂線の足からの距離の比が一定である曲面のことです。

楕円錐は$a=b=c$のとき球錐(単に「錐体」)、$a=b≠c$のとき回転楕円体、$a≠b≠c$のときは一般の楕円錐になります。

このように、楕円錐の方程式は比の関係で表され、その形状は$a$, $b$, $c$の値で決まります。

ソースコード

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
from mpl_toolkits.mplot3d import Axes3D

# パラメータ
a = 2  # x軸の半長径
b = 3  # y軸の半長径
c = 4  # z軸の半長径
u = 1  # x軸の平行移動量
v = 1  # y軸の平行移動量
w = 1  # z軸の平行移動量

# メッシュグリッドデータ
x = np.arange(-6, 6, 0.25)
y = np.arange(-6, 6, 0.25)
x, y = np.meshgrid(x, y)

z1 = (c / a) * np.sqrt(a**2 - (x - u)**2 - (y - v)**2) + w
z2 = -(c / a) * np.sqrt(a**2 - (x - u)**2 - (y - v)**2) + w

# 3Dプロット
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x, y, z1, alpha=0.5)
ax.plot_surface(x, y, z2, alpha=0.5)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()

このコードは、以下の手順で楕円錐のグラフを描画します。

1. NumPyとMatplotlibをインポートします。
2. 楕円錐のパラメータ(半長径と平行移動量)を設定します。
3. メッシュグリッドデータ($x$, $y$座標)を生成します。
4. 楕円錐の方程式から、$z$座標を計算します。
5. 3D空間にプロットし、軸ラベルを設定します。
6. グラフを表示します。

このコードを実行すると、指定された楕円錐の3Dグラフが表示されます。

パラメータを変更すれば、異なる楕円錐のグラフを描くことができます。

結果解説

[実行結果]

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

1. 表示内容

• 3次元の楕円錐の表面が、2つの半透明な曲面として描画されています。
• 1つの曲面は楕円錐の上半分を、もう1つの曲面は下半分を表しています。

2. 座標軸・ラベル

• 3次元の座標軸($X$軸、$Y$軸、$Z$軸)が表示されています。
• 各軸にはラベル(X、Y、Z)が付いています。

3. 曲面の形状

• 曲面の形状は、設定したパラメータ$(a, b, c, u, v, w)$によって決まります。
• $a$は$X$軸方向の半長径、$b$は$Y$軸方向の半長径、$c$はZ軸方向の半長径です。
• $u$、$v$、$w$は平行移動量で、楕円錐の頂点の位置を決めます。

4. 半透明効果

• 両方の曲面は$alpha=0.5$と設定されているため、半透明になっています。
• これにより、2つの曲面が重なった部分の形状がより分かりやすくなります。

5. 曲面の色

• デフォルトでは、曲面の色はオレンジ色とブルーの濃淡で表現されています。
• 色は曲面の高さ($z$座標の値)に応じて変化します。

このグラフは、与えられた楕円錐の方程式を視覚化したものです。

パラメータを変更すれば、異なる形状や位置の楕円錐を表すことができます。

半透明の重ね描画により、立体的な形状がよりわかりやすくなっています。

カドモフ方程式

カドモフ方程式（Kadomtsev-Petviashvili equation）は非線形偏微分方程式の一種であり、主に2次元非線形波動の理論において重要な方程式です。

この方程式は水理学やプラズマ物理学、非線形光学、または海洋学などさまざまな分野で応用されます。

一般的なカドモフ方程式は次の形を取ります：

$$
u_t + u_{xxx} + 6uu_x = 0
$$

ここで各項の意味は以下の通りです：

• $( u(x, t) ) $は変数$ ( x ) $と$ ( t ) $の関数であり、波の高さや系の状態を表します。
• $( u_t ) $は$ ( u ) $の時間微分であり、時間変化を表します。
• $( u_{xxx} ) $は$ ( u ) $の$ ( x ) $に関する3階微分であり、空間的な変化を表します。
• $( u_x ) $は$ ( u ) $の$ ( x ) $に関する1階微分であり、空間的な変化率を表します。
• $( u_{xxx} ) $は$ ( u_x ) $の$ ( x ) $に関する3階微分であり、空間的な変化率の空間微分を表します。
• $( u \cdot u_x ) $は非線形項であり、波の非線形相互作用を表します。

カドモフ方程式はソリトン（孤立波）やその他の非線形波動解を持つことで知られています。

この方程式は安定な非線形波動を記述し、系のエネルギー保存則や非線形な波動相互作用を理解するための基本的なモデルとして広く研究されています。

カドモフ方程式の解析的な解は一般に困難ですが、数値シミュレーションや解析的手法を用いて、方程式の振る舞いや特性を理解することが試みられています。

ソースコード

以下に、Pythonを使用してカドモフ方程式を解き、グラフ化する例を示します。

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

1
2
3
4

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

次に、カドモフ方程式の右辺を定義します。

1
2
3
4
5

def kadomtsev_petviashvili(t, u):
    u1, u2 = u
    du1_dt = u2
    du2_dt = 6 * u1 * (u1 - u2) * (1 - u1)
    return [du1_dt, du2_dt]

解析に使用する初期値と時間範囲を設定します。

1
2

t_span = [0, 10]
u0 = [0.5, 0.5]  # 初期値 [u1(0), u2(0)]

solve_ivp() を使用してカドモフ方程式を数値的に解きます。

1
2

sol = solve_ivp(kadomtsev_petviashvili, t_span, u0, t_eval=np.linspace(t_span[0], t_span[1], 1000))
u1, u2 = sol.y

最後に、解をグラフ化します。

1
2
3
4
5
6
7
8
9

plt.figure(figsize=(8, 6))
plt.plot(sol.t, u1, label='u1(t)')
plt.plot(sol.t, u2, label='u2(t)')
plt.xlabel('t')
plt.ylabel('u')
plt.title('Solution of the Kadomtsev-Petviashvili Equation')
plt.legend()
plt.grid(True)
plt.show()

このコードでは、solve_ivp() を使用してカドモフ方程式を解き、その解を時刻 t に対してプロットしています。

適切な初期値や時間範囲、グラフの設定を調整して、カドモフ方程式の数値解をグラフ化することができます。

[実行結果]

結果解説

カドモフ方程式の数値解をグラフ化する際、以下の内容が表示されます：

1. 軸のラベル:

• 横軸（x軸）: 時間 t
• 縦軸（y軸）: 解 u1(t) および u2(t) の値

2. 解の変化:

• グラフ上には時間 t に対する解 u1(t) と u2(t) の値がプロットされます。
• u1(t) と u2(t) は時間の関数として表現され、時間が増加するにつれてそれぞれの解の値がどのように変化するかが示されます。

3. 解の意味:

• u1(t) はカドモフ方程式の主要な解であり、系の状態や振る舞いを示します。
• u2(t) は u1(t) の時間微分に相当し、系の速度や変化率を表します。

4. タイトルと凡例:

• グラフには「Solution of the Kadomtsev-Petviashvili Equation」というタイトルが表示されます。
• 凡例には u1(t) と u2(t) のラベルが含まれ、どちらがどの曲線を表しているかが示されます。

5. グリッド:

• グラフには背景にグリッドが表示され、各目盛りが値の参照を補助します。

このグラフを通じて、時間とともにカドモフ方程式の解がどのように振る舞うかを視覚的に理解することができます。

u1(t) と u2(t) の値の変化を見ながら、方程式が記述する現象や振る舞いを把握することができます。

パラボロイドの方程式

パラボロイドは数学的な曲面であり、放物線を回転させて得られる回転放物面の一種です。
パラボロイドは二次曲面の一種であり、一般的には次のような形式で表されます：

$$
z = \frac{x^2}{a^2} + \frac{y^2}{b^2}
$$

ここで、$( x )$、$( y )$、$( z ) $は空間内の座標を表し、$( a ) $と$ ( b ) $は定数であり、それぞれ$ x軸$方向と$ y軸$方向の放物線の広がりを制御します。

この方程式における放物面は、$( x ) $と$ ( y ) $の値によって$ ( z ) $座標が決まり、放物線の頂点が原点に位置します。
パラボロイドは、放物線が$ x軸$と$ y軸$の両方向に広がっている形状を持ち、$ ( z ) $の値は放物線の広がりに応じて変化します。

Pythonでパラボロイドをグラフ化する際には、このような方程式を用いて$ ( x )$、$( y ) $の範囲を指定し、それに対応する$ ( z ) $の値を計算します。

得られた$ ( x )$、$( y )$、$( z ) $の値を3次元グラフとしてプロットすることで、パラボロイドの形状を視覚化することができます。

ソースコード

以下のPythonコードは、パラボロイドの方程式を解いてグラフ化するものです。

この例では、NumPyと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
36

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

# パラメータ設定
a = 1.0
b = 1.0
c = 1.0

# パラメータの範囲設定
u = np.linspace(-1, 1, 100)
v = np.linspace(-1, 1, 100)
U, V = np.meshgrid(u, v)

# パラボロイドの方程式をパラメータ表示で定義
X = a * U
Y = b * V
Z = c * (U**2 + V**2)

# 3Dプロット
fig = plt.figure(figsize=(8, 8))
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')

# グラフのタイトル
plt.title('Paraboloid')

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

このコードでは、a、b、cはパラメータとしてパラボロイドの形状を制御します。

UとVはパラメータの範囲を定義し、X、Y、Zはパラボロイドの方程式に基づいて表面の座標を計算します。

plot_surface関数を使用してパラボロイドを3Dプロットし、軸ラベルやタイトルを設定してグラフを表示します。

[実行結果]

ソースコード解説

このPythonコードは、パラボロイド曲面を生成し、3次元プロットで視覚化するためのものです。

以下にそれぞれの部分を詳しく説明します。

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

1
2
3

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

• numpyは数値計算を行うためのライブラリであり、配列や数学関数を効率的に扱います。
• matplotlib.pyplotはグラフ描画ライブラリで、グラフのプロットやカスタマイズに使用します。
• Axes3Dは3次元プロット用のクラスを提供します。

2. パラメータ設定:

1
2
3

a = 1.0
b = 1.0
c = 1.0

• パラボロイドの方程式における定数$ (a)$、$(b)$、$(c) $を設定します。

3. パラメータの範囲設定:

1
2
3

u = np.linspace(-1, 1, 100)
v = np.linspace(-1, 1, 100)
U, V = np.meshgrid(u, v)

• パラメータ$ (u)$、$(v) $の範囲を設定し、それらの値のメッシュグリッドを作成します。

4. パラボロイドの方程式を定義:

1
2
3

X = a * U
Y = b * V
Z = c * (U**2 + V**2)

• パラボロイドの方程式$ (z = c(u^2 + v^2)) $をパラメータ$ (U)$、$(V) $を使って定義します。

5. 3Dプロット:

1
2
3

fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, cmap='viridis', edgecolor='k')

• plot_surfaceメソッドを使用して、パラボロイド曲面を3次元プロットします。
• cmap='viridis'はカラーマップを設定し、’viridis’は色のマッピングを示します。
• edgecolor='k'は曲面の枠線の色を黒に設定します。

6. 軸ラベルとタイトルの設定:

1
2
3
4

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.title('Paraboloid')

• 軸ラベルとグラフのタイトルを設定します。

7. グラフの表示:

1

plt.show()

• 最後に、プロットされたグラフを表示します。

このコードを実行すると、$(u)$、$(v) $の範囲内で定義されたパラボロイド曲面が3次元空間に表示されます。

結果解説

[実行結果]

上記のPythonコードによって生成されるグラフは、パラボロイド（放物体）の3次元プロットです。

このグラフは以下のような特徴を持っています：

1. 形状:

• パラボロイドは放物線を回転軸周りに回転させて生成される曲面です。
  この例では、パラメータ a、b、c を調整することでパラボロイドの形状を変えることができます。

2. 色:

• cmap='viridis' を指定することで、グラフの色をViridisカラーマップに設定しています。
  これにより、高さ（Z軸方向）に応じて色が異なります。
  Viridisは高さが低いときは青色から、高いときは黄色や緑色に変化するカラーマップです。

3. 座標軸:

• グラフには$X軸$、$Y軸$、$Z軸$があります。
  各軸はそれぞれ X、Y、Z に対応し、パラボロイドの方程式の座標値を示しています。

4. タイトルとラベル:

• グラフには適切なタイトル ('Paraboloid') が付けられています。
• $X軸$、$Y軸$、$Z軸$にはそれぞれ 'X'、'Y'、'Z' というラベルが付けられており、各軸の意味を明示しています。

このグラフはパラメータを調整することでパラボロイドの形状を変えることができ、3次元空間で放物体の曲面を視覚化するのに役立ちます。

リャプノフ方程式は、非線形の偏微分方程式の一種で、以下のように表されます。

$$
∂u/∂t = -α u ∂u/∂x + β ∂^2u/∂x^2 + γ ∂(u^2)/∂x
$$

ここで、$u$ は未知関数、$t$ は時間、$x$ は空間座標、$α$、$β$、$γ$ はパラメータです。

この方程式は、様々な分野で現れる非線形現象をモデル化するのに役立ちます。
特に流体力学や プラズマ物理学、光学などの分野で重要な役割を果たします。

第1項　$ -α u ∂u/∂x $は移流項と呼ばれ、流れによる物理量の移動を表します。
第2項　$ β ∂^2u/∂x^2 $は拡散項で、物理量の空間的な広がりを表します。
第3項　$ γ ∂(u^2)/∂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

import numpy as np
import matplotlib.pyplot as plt

# パラメータの設定
alpha = 1.0
beta = 0.5
gamma = -1.0

# 時間と空間の範囲
tmin, tmax, dt = 0, 20, 0.01
xmin, xmax, dx = -10, 10, 0.1

# 初期条件
u0 = np.exp(-(x**2)/2)

# 時間と空間の格子点
t = np.arange(tmin, tmax, dt)
x = np.arange(xmin, xmax, dx)
nt, nx = len(t), len(x)

# 解の格納配列
u = np.zeros((nt, nx))

# 初期条件の設定
u[0] = u0

# Lyapunov 方程式の数値解法
for i in range(nt-1):
    for j in range(1, nx-1):
        u[i+1, j] = u[i, j] - alpha*dt*u[i, j]*(u[i, j+1] - u[i, j-1])/(2*dx) \
                    + beta*dt*(u[i, j+1] - 2*u[i, j] + u[i, j-1])/(dx**2) \
                    + gamma*dt*(u[i, j+1]**2 - u[i, j-1]**2)/(2*dx)

# グラフ描画
plt.figure(figsize=(10, 6))
plt.contourf(x, t, u, levels=np.linspace(np.min(u), np.max(u), 21), cmap='viridis')
plt.colorbar()
plt.xlabel('x')
plt.ylabel('t')
plt.title('Lyapunov Equation Solution')
plt.show()

上記のコードは、Lyapunov方程式の数値解を計算し、その結果をグラフ化するPythonプログラムです。

1. まず必要なライブラリ(numpy, matplotlib)をインポートします。
2. 次にLyapunov方程式のパラメータ(alpha, beta, gamma)を設定します。
3. 時間と空間の範囲を設定します。
4. 初期条件$u0$を定義します。
5. 時間と空間の格子点を作成します。
6. 解の格納配列$u$を初期化します。
7. 初期条件を設定します。
8. Lyapunov方程式の数値解法の本体部分です。
   時間ループと$x$方向の空間ループを使って、格子点上の解を計算していきます。
9. 最後に、matplotlibを使って計算結果をコンター図に描画します。

このプログラムを実行すると、Lyapunov方程式の数値解がコンター図として表示されます。
初期条件の影響を受けながら、時間の経過とともに解の様子が変化している様子が確認できます。

[実行結果]

ソースコード解説

ソースコードを詳しく説明します。

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

1
2

import numpy as np
import matplotlib.pyplot as plt

• Numpyライブラリと、Matplotlibライブラリをインポートしています。
• Numpyは数値計算に関する処理を行うのに使用します。
• Matplotlibは数値データの可視化を行うのに使用します。

2. パラメータの設定

1
2
3

alpha = 1.0
beta = 0.5 
gamma = -1.0

• リャプノフ方程式に現れるパラメータ$α$、$β$、$γ$の値を設定しています。

3. 時間と空間の範囲設定

1
2

tmin, tmax, dt = 0, 20, 0.01
xmin, xmax, dx = -10, 10, 0.1

• 時間$t$と空間$x$の計算範囲を設定しています。
• $t$は$0$から$20$まで、$dt=0.01$間隔で計算します。
• $x$は$-10$から$10$まで、$dx=0.1$間隔で計算します。

4. 初期条件の設定

1

u0 = np.exp(-(x**2)/2)

• 初期条件$u0$をガウス分布$exp(-(x**2)/2)$で与えています。

5. 時間と空間の格子点作成

1
2
3

t = np.arange(tmin, tmax, dt)
x = np.arange(xmin, xmax, dx)
nt, nx = len(t), len(x)

• 時間$t$と空間$x$の格子点を作成しています。
• 格子点の個数$nt$、$nx$も求めています。

6. 解の格納配列の初期化

1

u = np.zeros((nt, nx))

• 計算する解$u$を格納する配列を、$nt×nx$のゼロ行列として初期化しています。

7. 初期条件の代入

1

u[0] = u0 

• 解の格納配列の最初の時間ステップ$(t=0)$に初期条件$u0$を代入しています。

8. リャプノフ方程式の数値解法

1
2
3
4
5

for i in range(nt-1):
    for j in range(1, nx-1):
        u[i+1, j] = u[i, j] - alpha*dt*u[i, j]*(u[i, j+1] - u[i, j-1])/(2*dx) \
                    + beta*dt*(u[i, j+1] - 2*u[i, j] + u[i, j-1])/(dx**2) \
                    + gamma*dt*(u[i, j+1]**2 - u[i, j-1]**2)/(2*dx)

• 時間と空間の2重ループを使って、リャプノフ方程式の数値解法で次の時間ステップの解を計算しています。
• 移流項、拡散項、非線形項を陽的に計算しています。

9. グラフ描画

1
2
3
4
5
6
7

plt.figure(figsize=(10, 6))
plt.contourf(x, t, u, levels=np.linspace(np.min(u), np.max(u), 21), cmap='viridis')
plt.colorbar()
plt.xlabel('x')
plt.ylabel('t')
plt.title('Lyapunov Equation Solution')
plt.show()

• 最後に、Matplotlibを使って計算結果をコンター図で描画しています。
• 縦軸が時間$t$、横軸が空間$x$で、コンターの色分けが解$u$の値を表しています。
• 図のタイトル、軸ラベル、カラーバーも設定しています。

このように、このコードは初期条件から出発し、リャプノフ方程式の数値解法で次の時間ステップの解を計算し、最終的にその結果をコンター図で可視化するプログラムになっています。

結果解説

[実行結果]

このグラフは、時間$t$と空間$x$の2次元のコンター図となっています。
縦軸が$t$で時間の経過を表し、横軸が$x$で空間的な広がりを表しています。

コンター図の色分けは、方程式の解$u$の値を表しています。
コンターの色が赤に近づくほど解の値が大きく、青に近づくほど解の値が小さいことを意味しています。

時間$t=0$における初期条件は、$exp(-(x**2)/2)$によるガウス分布の形状になっていて、グラフの一番下の部分でそれが確認できます。
この初期条件から出発し、時間が経過するにつれて解の様子が変化していきます。

具体的には、初期条件のガウス分布が徐々に変形していき、時間ととともに特徴的なパターンが現れてきます。
このパターンの具体的な形状は、与えられた方程式のタイプや係数値によって異なります。

一般に、このようなコンター図を描くことで、与えられた非線形偏微分方程式の数値解がどのように時間と空間に対して進行するかを視覚的に捉えることができます。
初期条件の影響を受けながらも、独自のダイナミクスに従って解が変化する様子が確認できるのが特徴です。

このように、コンター図から方程式の解の時空間的な振る舞いを詳細に観察でき、対象となる現象の本質的な構造を理解する上で重要な情報を得ることができます。