既知の熱情報から分布を計算する

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

def solve_heat_distribution(known_temperatures, max_iter=1000, tolerance=1e-6,
                            grid_size=7):
    """
    5×5メッシュ上の熱分布を解く関数

    Parameters:
    -----------
    known_temperatures : dict
        既知の温度を持つセルの位置と温度値 {(i, j): temperature}
    max_iter : int
        最大反復回数
    tolerance : float
        収束判定の閾値

    Returns:
    --------
    temperature_grid : numpy.ndarray
        計算された温度分布
    iterations : int
        収束までの反復回数
    """
    # 5×5メッシュの初期化
    temperature_grid = np.zeros((grid_size, grid_size))

    # 既知の温度値を設定
    for (i, j), temp in known_temperatures.items():
        temperature_grid[i, j] = temp

    # 既知の位置のマスク（True=既知、False=未知）
    mask = np.zeros((grid_size, grid_size), dtype=bool)
    for i, j in known_temperatures.keys():
        mask[i, j] = True

    # 反復計算（ガウス・ザイデル法）
    iterations = 0
    for iter_count in range(max_iter):
        old_temperatures = temperature_grid.copy()

        # 各グリッドポイントを更新
        for i in range(grid_size):
            for j in range(grid_size):
                if not mask[i, j]:  # 未知の温度点のみ更新
                    # 隣接セルの温度を取得（境界を考慮）
                    neighbors = []
                    if i > 0:  # 上
                        neighbors.append(temperature_grid[i-1, j])
                    if i < grid_size-1:  # 下
                        neighbors.append(temperature_grid[i+1, j])
                    if j > 0:  # 左
                        neighbors.append(temperature_grid[i, j-1])
                    if j < grid_size-1:  # 右
                        neighbors.append(temperature_grid[i, j+1])

                    # 隣接セルの平均を計算
                    temperature_grid[i, j] = sum(neighbors) / len(neighbors)

        # 収束判定
        diff = np.max(np.abs(temperature_grid - old_temperatures))
        iterations = iter_count + 1

        if diff < tolerance:
            print(f"収束しました（反復回数: {iterations}）")
            break

    return temperature_grid, iterations

def visualize_temperature(temperature_grid, known_positions=None, grid_size=7):
    """
    温度分布を可視化する関数

    Parameters:
    -----------
    temperature_grid : numpy.ndarray
        計算された温度分布
    known_positions : list
        既知の温度を持つセルの位置 [(i, j), ...]
    """
    fig, ax = plt.subplots(figsize=(10, 8))

    # ヒートマップの描画
    im = ax.imshow(temperature_grid, cmap=cm.jet)

    # カラーバーの追加
    cbar = fig.colorbar(im, ax=ax)
    cbar.set_label('温度')

    # グリッド線の表示
    ax.set_xticks(np.arange(-.5, grid_size, 1), minor=True)
    ax.set_yticks(np.arange(-.5, grid_size, 1), minor=True)
    ax.grid(which="minor", color="black", linestyle='-', linewidth=2)

    # 座標表示
    ax.set_xticks(np.arange(0, grid_size, 1))
    ax.set_yticks(np.arange(0, grid_size, 1))

    # 既知の温度位置をマーク
    if known_positions:
        for i, j in known_positions:
            ax.text(j, i, f"{temperature_grid[i, j]:.1f}", ha="center", va="center",
                    color="black", fontweight="bold", bbox=dict(facecolor='white', alpha=0.7))

    # 値の表示
    for i in range(grid_size):
        for j in range(grid_size):
            if known_positions and (i, j) not in known_positions:
                ax.text(j, i, f"{temperature_grid[i, j]:.1f}", ha="center", va="center")

    plt.title('5×5メッシュの熱分布')
    plt.show()

# 使用例
if __name__ == "__main__":
    grid_size = 7

    # 既知の温度値を設定（例）
    known_temperatures = {
        (1, 1): 100,  # 左上角100度
        (2, 1): 80,  # 右上角: 100度
        (3, 1): 50,  # 右上角: 100度
        (4, 1): 80,  # 右上角: 100度
        (5, 1): 100,  # 右上角: 100度
        (1, 3): 30,    # 左下角: 0度
        (2, 3): 40,    # 右下角: 0度
        # (3, 2): 30,   # 中央: 50度
    }

    # 熱分布を計算
    temperature_grid, iterations = solve_heat_distribution(known_temperatures)

    # 結果を表示
    print("計算された温度分布:")
    print(temperature_grid)
    print(f"反復回数: {iterations}")

    # 結果を可視化
    visualize_temperature(temperature_grid, list(known_temperatures.keys()),
                          grid_size)

    # # 直接法（行列解法）による解法
    # print("\n直接法による解法:")

    # def solve_direct_method(known_temperatures, grid_size=5):
    #     """行列方程式を解いて熱分布を計算する関数"""
    #     # 未知の温度ポイントのインデックスを取得
    #     unknown_indices = []
    #     for i in range(grid_size):
    #         for j in range(grid_size):
    #             if (i, j) not in known_temperatures:
    #                 unknown_indices.append((i, j))

    #     n_unknowns = len(unknown_indices)

    #     # 係数行列と右辺ベクトルを構築
    #     A = np.zeros((n_unknowns, n_unknowns))
    #     b = np.zeros(n_unknowns)

    #     # インデックスマッピングを作成
    #     idx_map = {idx: k for k, idx in enumerate(unknown_indices)}

    #     # 方程式を構築
    #     for k, (i, j) in enumerate(unknown_indices):
    #         # 当該セルに対応する行に1を設定
    #         A[k, k] = 1

    #         # 隣接セルの貢献を計算
    #         neighbors = []
    #         neighbor_indices = []

    #         if i > 0:  # 上
    #             neighbors.append((i-1, j))
    #         if i < grid_size-1:  # 下
    #             neighbors.append((i+1, j))
    #         if j > 0:  # 左
    #             neighbors.append((i, j-1))
    #         if j < grid_size-1:  # 右
    #             neighbors.append((i, j+1))

    #         for ni, nj in neighbors:
    #             if (ni, nj) in known_temperatures:
    #                 # 既知の境界点からの貢献
    #                 b[k] += known_temperatures[(ni, nj)] / len(neighbors)
    #             else:
    #                 # 未知点からの貢献
    #                 nk = idx_map[(ni, nj)]
    #                 A[k, nk] = -1.0 / len(neighbors)

    #     # 連立方程式を解く
    #     x = np.linalg.solve(A, b)

    #     # 結果を元のグリッドに戻す
    #     result = np.zeros((grid_size, grid_size))
    #     for (i, j), temp in known_temperatures.items():
    #         result[i, j] = temp

    #     for k, (i, j) in enumerate(unknown_indices):
    #         result[i, j] = x[k]

    #     return result

    # # 直接法で解く
    # direct_solution = solve_direct_method(known_temperatures)
    # print(direct_solution)

    # # 直接法の結果を可視化
    # plt.figure(figsize=(10, 8))
    # plt.title('直接法による5×5メッシュの熱分布')
    # visualize_temperature(direct_solution, list(known_temperatures.keys()))