偏微分方程式(ラプラス方程式の離散化)

Pythonを用いて、与えられた温度分布を数値的に解き、ラプラス方程式の解を三次元プロットを作成するコード

Pythonコード

。

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

# 解析解を定義
def temperature_distribution(x, y):
    # 与えられた解析解 u(x, y) = (sin(πx) * sinh(π(1-y))) / sinh(π)
    pi = np.pi
    return (np.sin(pi * x) * np.sinh(pi * (1 - y))) / np.sinh(pi)

# メッシュグリッドを作成
x = np.linspace(0, 1, 100)
y = np.linspace(0, 1, 100)
X, Y = np.meshgrid(x, y)

# 温度分布を計算
Z = temperature_distribution(X, Y)

# プロットの設定
fig = plt.figure(figsize=(10, 8))
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('Temperature u(x, y)')
ax.set_title('Temperature Distribution on a Square Plate')

# カラーバーの追加
mappable = ax.plot_surface(X, Y, Z, cmap='viridis')
fig.colorbar(mappable, ax=ax, shrink=0.5, aspect=5)

# プロット表示
plt.show()

コードの説明

1. 温度分布の解析解:

   • temperature_distribution(x, y) 関数は、画像に示された解析解を用いて定義されています。
   • 温度分布は以下の式で与えられています：
     [
     u(x, y) = \frac{\sin(\pi x) \sinh(\pi (1-y))}{\sinh(\pi)}
     ]

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

   • x と y の範囲を 0 から 1 に分割し、100×100 のメッシュグリッドを作成します。

3. 温度分布の計算:

   • temperature_distribution(X, Y) を用いて、各メッシュポイントでの温度 Z を計算します。

4. 三次元プロット:

   • plot_surface 関数を使い、X, Y, Z を用いて三次元の温度分布をプロットします。
   • カラーマップ viridis を使用し、温度分布の視覚的な変化を示します。

5. カラーバーの追加:

   • プロットにカラーバーを追加し、温度分布の色と数値を対応させて視覚的に確認できるようにします。

結果の確認ポイント

• プロットでは、正方形板上の温度分布を三次元グラフとして確認できます。
• 温度分布は、周囲がゼロ温度に固定されており、境界条件を満たすように設定されています。

このコードを実行することで、板内部の温度分布を視覚化でき、ラプラス方程式の解析解がどのように板の内部温度を表現しているかを理解できます。必要に応じて、メッシュサイズやカラーマップの設定を変更し、異なる視覚表現を試すこともできます。

import numpy as np
import matplotlib.pyplot as plt

# 関数とその真の微分の定義
def f(x):
    return np.sin(x)  # 例として sin(x) を使用

def f_prime_exact(x):
    return np.cos(x)  # sin(x) の真の微分

# 差分近似の定義
def forward_difference(f, x, h):
    return (f(x + h) - f(x)) / h  # 前進差分

def backward_difference(f, x, h):
    return (f(x) - f(x - h)) / h  # 後退差分

def central_difference(f, x, h):
    return (f(x + h) - f(x - h)) / (2 * h)  # 中心差分

# 計算範囲とステップサイズの設定
x = np.linspace(0, 2 * np.pi, 100)  # 0 から 2π までの 100 ステップ
h = 0.1  # ステップサイズ

# 各差分法による微分近似
forward_diff = forward_difference(f, x, h)
backward_diff = backward_difference(f, x, h)
central_diff = central_difference(f, x, h)

# 真の微分の計算
exact_diff = f_prime_exact(x)

# 結果のプロット
plt.figure(figsize=(12, 8))

# 真の微分
plt.plot(x, exact_diff, label='True Derivative (cos(x))', linestyle='-', linewidth=2, color='black')

# 前進差分
plt.plot(x, forward_diff, label='Forward Difference Approximation', linestyle='--', color='blue')

# 後退差分
plt.plot(x, backward_diff, label='Backward Difference Approximation', linestyle='-.', color='green')

# 中心差分
plt.plot(x, central_diff, label='Central Difference Approximation', linestyle=':', color='red')

# グラフの設定
plt.xlabel('x')
plt.ylabel("f'(x)")
plt.title('Numerical Differentiation Using Difference Approximations')
plt.legend()
plt.grid()
plt.show()

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

# 2次元関数の定義
def u(x, y):
    return np.sin(np.pi * x) * np.sin(np.pi * y)

# 偏微分の差分近似を定義
def forward_difference_x(u, x, y, hx):
    return (u(x + hx, y) - u(x, y)) / hx

def backward_difference_x(u, x, y, hx):
    return (u(x, y) - u(x - hx, y)) / hx

def central_difference_x(u, x, y, hx):
    return (u(x + hx, y) - u(x - hx, y)) / (2 * hx)

def central_difference_y(u, x, y, hy):
    return (u(x, y + hy) - u(x, y - hy)) / (2 * hy)

def second_derivative_x(u, x, y, hx):
    return (u(x + hx, y) - 2 * u(x, y) + u(x - hx, y)) / hx**2

def second_derivative_y(u, x, y, hy):
    return (u(x, y + hy) - 2 * u(x, y) + u(x, y - hy)) / hy**2

# グリッドの設定
hx = 0.1  # x方向のステップサイズ
hy = 0.1  # y方向のステップサイズ
x = np.arange(0, 1 + hx, hx)
y = np.arange(0, 1 + hy, hy)
X, Y = np.meshgrid(x, y)

# 各点での関数値を計算
U = u(X, Y)

# 各差分法を用いて偏微分を計算
dU_dx_forward = forward_difference_x(u, X, Y, hx)
dU_dx_backward = backward_difference_x(u, X, Y, hx)
dU_dx_central = central_difference_x(u, X, Y, hx)
dU_dy_central = central_difference_y(u, X, Y, hy)

# 二階偏微分を中心差分で計算
d2U_dx2 = second_derivative_x(u, X, Y, hx)
d2U_dy2 = second_derivative_y(u, X, Y, hy)

# プロット設定
fig = plt.figure(figsize=(15, 10))

# 前進差分による1階偏微分
ax1 = fig.add_subplot(231, projection='3d')
ax1.plot_surface(X, Y, dU_dx_forward, cmap='viridis')
ax1.set_title("Forward Difference ∂u/∂x")

# 後退差分による1階偏微分
ax2 = fig.add_subplot(232, projection='3d')
ax2.plot_surface(X, Y, dU_dx_backward, cmap='plasma')
ax2.set_title("Backward Difference ∂u/∂x")

# 中心差分による1階偏微分 (x方向)
ax3 = fig.add_subplot(233, projection='3d')
ax3.plot_surface(X, Y, dU_dx_central, cmap='coolwarm')
ax3.set_title("Central Difference ∂u/∂x")

# 中心差分による1階偏微分 (y方向)
ax4 = fig.add_subplot(234, projection='3d')
ax4.plot_surface(X, Y, dU_dy_central, cmap='viridis')
ax4.set_title("Central Difference ∂u/∂y")

# 二階偏微分 (x方向)
ax5 = fig.add_subplot(235, projection='3d')
ax5.plot_surface(X, Y, d2U_dx2, cmap='plasma')
ax5.set_title("Second Derivative ∂²u/∂x²")

# 二階偏微分 (y方向)
ax6 = fig.add_subplot(236, projection='3d')
ax6.plot_surface(X, Y, d2U_dy2, cmap='coolwarm')
ax6.set_title("Second Derivative ∂²u/∂y²")

# プロット表示
plt.tight_layout()
plt.show()

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

# 領域の設定
N = 20  # グリッドの分割数 (内部点の数)
h = 1 / N  # グリッドの刻み幅
x = np.linspace(0, 1, N + 1)  # x軸方向の離散点
y = np.linspace(0, 1, N + 1)  # y軸方向の離散点
X, Y = np.meshgrid(x, y)

# 境界条件の設定
u = np.zeros((N + 1, N + 1))  # すべての点での温度分布の初期化
u[:, 0] = np.sin(np.pi * y)  # 左側境界 (u(x, 0) = sin(πy))
u[:, -1] = np.sin(np.pi * y)  # 右側境界 (u(x, 1) = sin(πy))
u[0, :] = 0  # 下側境界 (u(0, y) = 0)
u[-1, :] = 0  # 上側境界 (u(1, y) = 0)

# 内部点の反復法による解の計算 (ガウス・ザイデル法)
tolerance = 1e-5  # 収束条件
max_iterations = 5000  # 最大反復回数

for iteration in range(max_iterations):
    u_old = u.copy()
    
    # 内部点の更新
    for i in range(1, N):
        for j in range(1, N):
            u[i, j] = 0.25 * (u_old[i+1, j] + u_old[i-1, j] + u_old[i, j+1] + u_old[i, j-1])
    
    # 収束判定
    error = np.linalg.norm(u - u_old, ord=2)
    if error < tolerance:
        print(f"Converged after {iteration} iterations with error {error:.5e}.")
        break

# 結果のプロット
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')

# 温度分布の三次元プロット
ax.plot_surface(X, Y, u, cmap='viridis')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('Temperature u(x, y)')
ax.set_title('Solution of Laplace Equation Using Finite Difference Method')

# カラーバーの追加
mappable = ax.plot_surface(X, Y, u, cmap='viridis')
fig.colorbar(mappable, ax=ax, shrink=0.5, aspect=5)

# プロット表示
plt.show()