Pythonとコピペで学べる数値偏微分方程式

# プログラム名: fourier_2d_series_plot.py
# 内容: 2次元フーリエ級数を有限項まで展開して f(x, y) をプロット

import numpy as np
import matplotlib.pyplot as plt

L = 1  # 区間長
N = 20  # 項数
x = np.linspace(0, L, 100)
y = np.linspace(0, L, 100)
X, Y = np.meshgrid(x, y)

f_xy = np.zeros_like(X)

for m in range(1, N+1):
    for n in range(1, N+1):
        lambda_mn = (1 - np.cos(m * np.pi / 2)) * (1 - np.cos(n * np.pi / 2)) / (m * n)
        term = lambda_mn * np.sin(m * np.pi * X / L) * np.sin(n * np.pi * Y / L)
        f_xy += term

f_xy *= 4 / np.pi**2

# --- 可視化 ---
plt.figure(figsize=(6, 5))
cp = plt.contourf(X, Y, f_xy, levels=100, cmap='viridis')
plt.colorbar(cp)
plt.title("2D Fourier Series Approximation of f(x, y)")
plt.xlabel("x")
plt.ylabel("y")
plt.tight_layout()
plt.show()

# プログラム名: fourier_transform_table_plot_fixed_final.py
# 内容: 代表的な関数（矩形を除く）のフーリエ変換を可視化（エラー修正済み）

import numpy as np
import matplotlib.pyplot as plt
from sympy import symbols, exp, Abs, fourier_transform, Heaviside, DiracDelta, pi, sqrt, lambdify, simplify
from sympy.abc import x, alpha, p

# --- パラメータ設定 ---
p_val = 1

# --- 関数定義とフーリエ変換（解析） ---
f2 = exp(-p * Abs(x))                  
F2 = simplify(fourier_transform(f2, x, alpha))

f3 = DiracDelta(x)                     
F3 = simplify(fourier_transform(f3, x, alpha))

f4 = Heaviside(x)                      
F4 = simplify(fourier_transform(f4, x, alpha))

f5 = exp(-p * x**2)                    
F5 = simplify(fourier_transform(f5, x, alpha))

# --- lambdify で数値化（simplify後に限定） ---
a_vals = np.linspace(-10, 10, 400)
f2_lamb = lambdify(alpha, F2.subs(p, p_val), "numpy")
f3_lamb = lambdify(alpha, F3, "numpy")
f4_lamb = lambdify(alpha, F4, "numpy")
f5_lamb = lambdify(alpha, F5.subs(p, p_val), "numpy")

# --- グラフ描画 ---
plt.figure(figsize=(12, 8))

plt.subplot(2, 2, 1)
plt.plot(a_vals, f2_lamb(a_vals))
plt.title("Fourier Transform of Exponential Decay")
plt.xlabel("α")
plt.ylabel("F(α)")

plt.subplot(2, 2, 2)
plt.plot(a_vals, f3_lamb(a_vals))
plt.title("Fourier Transform of Delta Function")
plt.xlabel("α")
plt.ylabel("F(α)")

plt.subplot(2, 2, 3)
plt.plot(a_vals, np.real(f4_lamb(a_vals)))
plt.title("Fourier Transform of Heaviside Step")
plt.xlabel("α")
plt.ylabel("Re(F(α))")

plt.subplot(2, 2, 4)
plt.plot(a_vals, f5_lamb(a_vals))
plt.title("Fourier Transform of Gaussian Function")
plt.xlabel("α")
plt.ylabel("F(α)")

plt.tight_layout()
plt.show()

# プログラム名: fft_transform_table_plot.py
# 内容: 表の各 f(x) に対して離散フーリエ変換 (FFT) を実施し、スペクトルを可視化する

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq, fftshift

# --- 共通設定 / Common settings ---
N = 1024                   # サンプル数 / Number of samples
L = 10                     # 区間の長さ / Interval length
dx = L / N                 # サンプル間隔 / Sampling interval
x = np.linspace(-L/2, L/2, N, endpoint=False)
freqs = fftshift(fftfreq(N, d=dx))

# --- 1. 矩形関数（width=2r） / Rectangular pulse ---
r = 1
f1 = np.where(np.abs(x) <= r, 1/(2*r), 0)

# --- 2. 指数関数 / Exponential decay ---
p = 1
f2 = np.exp(-p * np.abs(x))

# --- 3. デルタ関数（近似: narrow peak） / Delta function approximation ---
f3 = np.zeros_like(x)
f3[N//2] = 1 / dx  # 中心にピーク

# --- 4. ヘヴィサイド関数 / Heaviside step function ---
f4 = np.heaviside(x, 1)

# --- 5. ガウス関数 / Gaussian ---
f5 = np.exp(-p * x**2)

# --- FFTを実行 / Compute FFT ---
def compute_fft(f):
    F = fft(f)
    F_shifted = fftshift(F)
    return np.abs(F_shifted)

F1 = compute_fft(f1)
F2 = compute_fft(f2)
F3 = compute_fft(f3)
F4 = compute_fft(f4)
F5 = compute_fft(f5)

# --- プロット / Plot results ---
plt.figure(figsize=(12, 10))

plt.subplot(3, 2, 1)
plt.plot(freqs, F1)
plt.title("FFT of Rectangular Pulse")

plt.subplot(3, 2, 2)
plt.plot(freqs, F2)
plt.title("FFT of Exponential Decay")

plt.subplot(3, 2, 3)
plt.plot(freqs, F3)
plt.title("FFT of Delta Function (Approx)")

plt.subplot(3, 2, 4)
plt.plot(freqs, F4)
plt.title("FFT of Heaviside Step")

plt.subplot(3, 2, 5)
plt.plot(freqs, F5)
plt.title("FFT of Gaussian Function")

plt.tight_layout()
plt.show()

# プログラム名: convolution_integral_numpy.py
# 内容: NumPyによる畳み込み積分（台形則で近似）

import numpy as np
import matplotlib.pyplot as plt

# --- 定義域と関数定義 / Define time domain and functions ---
t = np.linspace(0, 10, 1000)
f = np.exp(-t)          # f(t) = e^(-t)
g = np.heaviside(t, 1)  # g(t) = u(t)

# --- 畳み込み積分 / Convolution integral using trapezoidal rule ---
dt = t[1] - t[0]
conv = np.zeros_like(t)

for i in range(len(t)):
    integrand = f[:i+1] * g[i::-1]
    conv[i] = np.trapz(integrand, dx=dt)

# --- プロット / Plot result ---
plt.figure(figsize=(8, 4))
plt.plot(t, conv, label="f * g (convolution)")
plt.title("Convolution Integral: f(t) = e^(-t), g(t) = u(t)")
plt.xlabel("t")
plt.ylabel("(f * g)(t)")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
#wave_equation_fourier_solution.py


# --- パラメータ設定 ---
L = 1.0      # 弦の長さ / String length
v = 1.0      # 波の速度 / Wave speed
h = 1.0      # 振幅スケール
N = 50       # 展開項数 / Number of terms in series

# --- 座標と時間軸設定 ---
x = np.linspace(0, L, 500)
t_vals = [0, 0.2, 0.5, 1.0]  # 描画する時刻

# --- u(x, t) の定義 ---
def u_xt(x, t, L, v, h, N):
    result = np.zeros_like(x)
    for m in range(1, N+1):
        coef = (1 / m**2) * np.sin(m * np.pi / 3)
        result += coef * np.sin(m * np.pi * x / L) * np.cos(m * np.pi * v * t / L)
    return (9 * h / np.pi**2) * result

# --- グラフ描画 ---
plt.figure(figsize=(10, 6))
for t in t_vals:
    u = u_xt(x, t, L, v, h, N)
    plt.plot(x, u, label=f"t = {t:.2f}")

plt.title("Wave Equation Solution via Fourier Series")
plt.xlabel("x")
plt.ylabel("u(x, t)")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

# プログラム名: wave_2d_fourier_series.py
# 内容: 2次元波動方程式のフーリエ級数解を数値的に計算・プロットする
# Purpose: Visualize Fourier series solution to 2D wave equation on a unit square

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータ定義 / Define constants ---
v = 1          # 波の速度 / Wave speed
N = 20         # 項数 / Number of terms in the series

# --- 座標格子 / Spatial grid ---
x = np.linspace(0, 1, 100)
y = np.linspace(0, 1, 100)
X, Y = np.meshgrid(x, y)

# --- 計算対象の時間 / Target time
t = 0.3

# --- フーリエ級数による解 / Fourier series solution ---
def u_xy_t(X, Y, t, v, N):
    result = np.zeros_like(X)
    for m in range(1, N + 1):
        for n in range(1, N + 1):
            # 偶数項は打ち消される / Even terms cancel due to (1 - (-1)^m) = 0
            if m % 2 == 0 or n % 2 == 0:
                continue
            coef = (1 - (-1)**m) * (1 - (-1)**n) / (m**3 * n**3)
            result += coef * np.sin(m * np.pi * X) * np.sin(n * np.pi * Y) * \
                      np.cos(v * np.pi * np.sqrt(m**2 + n**2) * t)
    return (16 / np.pi**6) * result

# --- 計算 / Compute solution ---
U = u_xy_t(X, Y, t, v, N)

# --- 可視化 / Plot the solution ---
plt.figure(figsize=(6, 5))
cp = plt.contourf(X, Y, U, levels=50, cmap='plasma')
plt.colorbar(cp)
plt.title(f"2D Wave Equation u(x, y, t={t})")
plt.xlabel("x")
plt.ylabel("y")
plt.tight_layout()
plt.show()

# プログラム名: wave_2d_solution_given_form.py
# 内容: 与えられた 2D 波動方程式のフーリエ級数解を可視化
# Purpose: Visualize specific Fourier solution of the 2D wave equation

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータ設定 / Parameters ---
v = 1.0         # 波速 / wave speed
N = 20          # 項数 / number of terms
t = 0.5         # 時間 / time

# --- 座標グリッド / Spatial grid ---
x = np.linspace(0, 1, 100)
y = np.linspace(0, 1, 100)
X, Y = np.meshgrid(x, y)

# --- 解の定義 / Define the solution ---
def u_solution(X, Y, t, v, N):
    result = np.zeros_like(X)
    for m in range(1, N+1):
        for n in range(1, N+1):
            sin_mp2 = np.sin(m * np.pi / 2)
            sin_np2 = np.sin(n * np.pi / 2)
            term = (1 / (m**2 * n**2)) * sin_mp2 * sin_np2
            result += term * np.sin(m * np.pi * X / 2) * np.sin(n * np.pi * Y / 2) * \
                      np.cos((np.pi * v * t / 2) * np.sqrt(m**2 + n**2))
    return (64 / np.pi**4) * result

# --- 計算と可視化 / Compute and plot ---
U = u_solution(X, Y, t, v, N)

plt.figure(figsize=(6, 5))
cp = plt.contourf(X, Y, U, levels=50, cmap='viridis')
plt.colorbar(cp)
plt.title(f"2D Wave Solution at t = {t}")
plt.xlabel("x")
plt.ylabel("y")
plt.tight_layout()
plt.show()

# プログラム名: heat_eq_solution_fourier.py
# 内容: 一次元熱伝導方程式の解析解（フーリエ級数）を計算・可視化
# Purpose: Visualize Fourier series solution of 1D heat equation

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータ設定 / Parameters ---
L = 1.0         # 長さ / Length of rod
a = 0.1         # 熱拡散率 / Thermal diffusivity
T0 = 100        # 初期温度スケール
N = 50          # フーリエ項数 / Number of terms in series

# --- 空間・時間グリッド / Define spatial and temporal grids ---
x = np.linspace(0, L, 500)
t_vals = [0.01, 0.1, 0.5, 1.0]

# --- 解の定義 / Define the Fourier series solution ---
def u_xt(x, t, L, a, T0, N):
    result = np.zeros_like(x)
    for m in range(1, N + 1):
        coef = (1 / m) * (np.cos(m * np.pi / 2) - np.cos(3 * m * np.pi / 4))
        result += coef * np.sin(m * np.pi * x / L) * np.exp(- (m * np.pi / L)**2 * a * t)
    return (2 * T0 / np.pi) * result

# --- プロット / Plot the solution ---
plt.figure(figsize=(10, 6))
for t in t_vals:
    u = u_xt(x, t, L, a, T0, N)
    plt.plot(x, u, label=f"t = {t}")

plt.title("1D Heat Equation Solution via Fourier Series")
plt.xlabel("x")
plt.ylabel("u(x, t)")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# プログラム名: heat_2d_fourier_coefficients.py
# 内容: 2D熱伝導方程式の初期係数 A_mn を可視化
# Purpose: Visualize Fourier coefficients for initial 2D temperature profile

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータ設定 / Parameters ---
T0 = 100
M, N = 30, 30  # m, n の最大値

A = np.zeros((M, N))

# --- 係数の計算 / Compute Fourier coefficients A_mn ---
for m in range(1, M+1):
    for n in range(1, N+1):
        term_m = np.cos(m * np.pi / 2) - (-1)**m
        term_n = np.cos(n * np.pi / 2) - (-1)**n
        A[m-1, n-1] = (1 / (m * n)) * term_m * term_n

A *= (4 * T0) / (np.pi**2)

# --- 可視化 / Visualization ---
plt.figure(figsize=(6, 5))
plt.imshow(A, extent=[1, N, 1, M], origin='lower', aspect='auto', cmap='plasma')
plt.colorbar(label="Coefficient A[m,n]")
plt.title("Fourier Coefficients of Initial Condition")
plt.xlabel("n")
plt.ylabel("m")
plt.tight_layout()
plt.show()

# プログラム名: heat_eq_gaussian_solution.py
# 内容: 1次元熱方程式の基本解（グリーン関数）を描画
# Purpose: Plot the fundamental Gaussian solution to the 1D heat equation

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータ設定 / Parameters ---
a = 0.5                          # 熱拡散率 / Thermal diffusivity
x = np.linspace(-10, 10, 500)   # 空間範囲 / Spatial domain
t_vals = [0.01, 0.1, 0.5, 1.0]  # 表示する時間ステップ / Time steps

# --- 解の定義 / Define the fundamental solution u(x, t) ---
def u_xt(x, t, a):
    return (1 / (2 * np.sqrt(np.pi * a * t))) * np.exp(-x**2 / (4 * a * t))

# --- プロット / Plot ---
plt.figure(figsize=(10, 6))
for t in t_vals:
    plt.plot(x, u_xt(x, t, a), label=f"t = {t}")

plt.title("Fundamental Solution to 1D Heat Equation")
plt.xlabel("x")
plt.ylabel("u(x, t)")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

# プログラム名: heat_equation_nonhomogeneous.py
# 内容: 非斉次境界条件を持つ熱伝導方程式の解析解を描画
# Purpose: Visualize 1D heat equation with non-homogeneous boundary condition

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータ設定 / Parameters ---
L = 1.0         # 棒の長さ / Length of rod
a = 0.1         # 熱拡散係数 / Thermal diffusivity
T0 = 100        # 境界温度スケール / Boundary temperature
N = 50          # フーリエ展開項数 / Number of terms
x = np.linspace(0, L, 500)
t_vals = [0.01, 0.1, 0.5, 1.0]  # 時刻 / Time snapshots

# --- 解の定義 / Solution function u(x, t) ---
def u_xt(x, t, L, a, T0, N):
    result = np.zeros_like(x)
    for m in range(1, N+1):
        coef = (1 / m) * (-1)**m * (np.cos(3 * m * np.pi / 4) - np.cos(m * np.pi / 2))
        result += coef * np.sin(m * np.pi * x / L) * np.exp(- (m * np.pi / L)**2 * a * t)
    return (T0 / np.pi) * result + (T0 / (2 * L)) * x

# --- プロット / Plot ---
plt.figure(figsize=(10, 6))
for t in t_vals:
    u = u_xt(x, t, L, a, T0, N)
    plt.plot(x, u, label=f"t = {t}")

plt.title("1D Heat Equation with Nonhomogeneous Boundary Condition")
plt.xlabel("x")
plt.ylabel("u(x, t)")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

# プログラム名: laplace_2d_fourier_solution.py
# 内容: 2次元ラプラス方程式の解析解（フーリエ展開）を可視化
# Purpose: Visualize 2D Laplace equation solution via Fourier series

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータ設定 / Parameters ---
T0 = 100         # 左側境界温度 / Boundary temperature at x = 0
L = 1.0          # 領域サイズ / Domain length
N = 50           # フーリエ項数 / Number of terms

# --- グリッド / Grid setup ---
x = np.linspace(0, L, 100)
y = np.linspace(0, L, 100)
X, Y = np.meshgrid(x, y)

# --- 解の定義 / Solution function ---
def u_xy(X, Y, T0, L, N):
    result = np.zeros_like(X)
    for m in range(1, N+1):
        numerator = np.cos(m * np.pi / 2) - 1
        denom = m * np.sinh(m * np.pi)
        term = (numerator / denom) * np.sinh(m * np.pi * (X / L - 1)) * np.sin(m * np.pi * Y / L)
        result += term
    return (2 * T0 / np.pi) * result

# --- 計算とプロット / Compute and plot ---
U = u_xy(X, Y, T0, L, N)

plt.figure(figsize=(6, 5))
cp = plt.contourf(X, Y, U, levels=50, cmap='plasma')
plt.colorbar(cp)
plt.title("Laplace Equation Solution u(x, y)")
plt.xlabel("x")
plt.ylabel("y")
plt.tight_layout()
plt.show()

# プログラム名: laplace_solution_asymmetric.py
# 内容: 与えられた2Dラプラス方程式の特解を数値的に可視化
# Purpose: Visualize specific solution to 2D Laplace equation with asymmetric BC

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータ設定 / Parameters ---
T0 = 100
N = 50  # フーリエ項数 / Number of Fourier terms

# --- グリッド設定 / Grid for (x, y) ---
x = np.linspace(0, 1, 100)
y = np.linspace(0, 1, 100)
X, Y = np.meshgrid(x, y)

# --- 解の定義 / Define solution u(x, y) ---
def u_xy(X, Y, T0, N):
    result = np.zeros_like(X)
    for m in range(1, N + 1):
        if (m % 2) == 0:  # (-1)^m - 1 = 0 for even m
            continue
        coef = ((-1)**m - 1) / m  # This will be -2/m for odd m
        result += coef * (1 / np.sinh(m * np.pi / 2)) * \
                  np.sinh(m * np.pi * (X - 0.5)) * np.sin(m * np.pi * Y)
    return (2 * T0 / np.pi) * result

# --- 可視化 / Visualization ---
U = u_xy(X, Y, T0, N)

plt.figure(figsize=(6, 5))
cp = plt.contourf(X, Y, U, levels=50, cmap='viridis')
plt.colorbar(cp)
plt.title("Solution to Laplace Equation u(x, y)")
plt.xlabel("x")
plt.ylabel("y")
plt.tight_layout()
plt.show()

# プログラム名: bessel_j1_plot.py
# 内容: 第1種ベッセル関数 J_1(x) の数値計算とグラフ表示
# Purpose: Plot the first kind Bessel function J1(x) using scipy

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import jn  # jn(n, x) は Bessel 関数 J_n(x)

# --- x軸設定 / Define x range ---
x = np.linspace(0, 20, 500)

# --- J_1(x) の計算 / Compute J1(x) ---
J1 = jn(1, x)

# --- プロット / Plot ---
plt.figure(figsize=(8, 5))
plt.plot(x, J1, label=r"$J_1(x)$", color='blue')
plt.title("Bessel Function of the First Kind: $J_1(x)$")
plt.xlabel("x")
plt.ylabel("J₁(x)")
plt.grid(True)
plt.axhline(0, color='black', linewidth=0.8)
plt.legend()
plt.tight_layout()
plt.show()

# プログラム名: pde_classification_examples.py
# 内容: 熱方程式（放物型）、波動方程式（双曲型）、ラプラス方程式（楕円型）の解を可視化
# Purpose: Visual demonstration of PDE classification (Parabolic, Hyperbolic, Elliptic)

import numpy as np
import matplotlib.pyplot as plt

# --- 空間と時間のメッシュ / Create mesh grid ---
x = np.linspace(0, 1, 100)
t = np.linspace(0, 1, 100)
X, T = np.meshgrid(x, t)

# --- 放物型: 熱方程式の一例解（Fourier解） ---
U_heat = np.exp(-np.pi**2 * T) * np.sin(np.pi * X)

# --- 双曲型: 波動方程式の一例解（Standing wave） ---
U_wave = np.sin(np.pi * X) * np.cos(np.pi * T)

# --- 楕円型: ラプラス方程式の一例解（定常温度分布） ---
Y = np.linspace(0, 1, 100)
X_ell, Y_ell = np.meshgrid(x, Y)
U_laplace = np.sinh(np.pi * X_ell) * np.sin(np.pi * Y_ell)

# --- 可視化 / Plotting ---
fig, axs = plt.subplots(1, 3, figsize=(18, 5))

# Heat Equation (Parabolic)
cs1 = axs[0].contourf(X, T, U_heat, cmap='plasma')
axs[0].set_title("Parabolic: Heat Equation")
axs[0].set_xlabel("x")
axs[0].set_ylabel("t")
fig.colorbar(cs1, ax=axs[0])

# Wave Equation (Hyperbolic)
cs2 = axs[1].contourf(X, T, U_wave, cmap='viridis')
axs[1].set_title("Hyperbolic: Wave Equation")
axs[1].set_xlabel("x")
axs[1].set_ylabel("t")
fig.colorbar(cs2, ax=axs[1])

# Laplace Equation (Elliptic)
cs3 = axs[2].contourf(X_ell, Y_ell, U_laplace, cmap='coolwarm')
axs[2].set_title("Elliptic: Laplace Equation")
axs[2].set_xlabel("x")
axs[2].set_ylabel("y")
fig.colorbar(cs3, ax=axs[2])

plt.tight_layout()
plt.show()

# プログラム名: pde_analytical_methods_demo.py
# 内容: PDE解析解法の代表例（分離変数法、フーリエ展開、ラプラス変換的解釈、グリーン関数的解）を可視化
# Purpose: Demonstrate key PDE analytical techniques in one program

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import jn  # ベッセル関数 Jn(x)

# --- 共通定義 / Common domain ---
x = np.linspace(0, 1, 500)
t = np.linspace(0, 1, 500)
X, T = np.meshgrid(x, t)

# === 1. 分離変数法 + フーリエ級数による熱方程式の解 ===
U_heat = np.exp(-np.pi**2 * T) * np.sin(np.pi * X)

# === 2. ラプラス変換的な視点での指数減衰応答 ===
U_laplace = (1 - np.exp(-5 * T)) * np.heaviside(T, 1)  # ステップ応答的な振る舞い

# === 3. グリーン関数解（ガウス関数）===
def green_func(x, t, a=0.1):
    return (1 / (2 * np.sqrt(np.pi * a * t))) * np.exp(-x**2 / (4 * a * t))
G = green_func(X - 0.5, T + 0.01)  # シフト＆微小時間から開始

# === 4. 特異関数展開（ベッセル関数 J1(x)） ===
x_bessel = np.linspace(0, 20, 500)
J1 = jn(1, x_bessel)

# --- 可視化 / Visualization ---
fig, axs = plt.subplots(2, 2, figsize=(14, 10))

# 1. 熱方程式：分離変数 + フーリエ級数
axs[0, 0].contourf(X, T, U_heat, cmap='plasma')
axs[0, 0].set_title("1. Heat Equation (Separation + Fourier)")
axs[0, 0].set_xlabel("x")
axs[0, 0].set_ylabel("t")

# 2. ラプラス変換視点
axs[0, 1].plot(t, U_laplace[250], label="Step Response")
axs[0, 1].set_title("2. Laplace-like Response")
axs[0, 1].set_xlabel("t")
axs[0, 1].set_ylabel("u(t)")
axs[0, 1].grid(True)

# 3. グリーン関数（ガウス解）
axs[1, 0].contourf(X, T, G, cmap='viridis')
axs[1, 0].set_title("3. Green's Function: Heat Kernel")
axs[1, 0].set_xlabel("x")
axs[1, 0].set_ylabel("t")

# 4. ベッセル関数 J1
axs[1, 1].plot(x_bessel, J1, color='blue')
axs[1, 1].set_title("4. Bessel Function $J_1(x)$")
axs[1, 1].set_xlabel("x")
axs[1, 1].set_ylabel("J₁(x)")
axs[1, 1].grid(True)

plt.tight_layout()
plt.show()

# プログラム名: pde_analytical_methods_combined.py
# 内容: PDE解析解法・変換手法（分離変数法＋フーリエ級数、グリーン関数、フーリエ変換、特異関数展開）の例を統合し、
#       数値計算と可視化を行う
# Purpose: Demonstrate key analytical methods for PDEs including separation of variables with Fourier series,
#          Green's function, Fourier transform method, and special function expansion (Bessel functions).

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import jn  # Bessel function of the first kind
from scipy.fft import fft, fftshift, fftfreq

### 1. Separation of Variables & Fourier Series for Heat Equation
# ---------------------------------------------------------------
# PDE: ∂u/∂t = a ∂²u/∂x²,  0 < x < L, t > 0, with u(0,t)=u(L,t)=0
# Initial condition u(x,0) = sin(π x)  (for simplicity)
# Analytical solution is approximated by:
#   u(x,t) = sum_{n=1}^{N} b_n * sin(nπx/L) * exp(-(nπ/L)² a t)
#
# In this example, with initial condition sin(π x), only n=1 coefficient is nonzero.
a = 0.01          # thermal diffusivity
L = 1.0           # length of domain
Nx = 200
x = np.linspace(0, L, Nx)
t_heat = 0.5      # evaluation time

# Fourier series solution (only fundamental mode active)
N_modes = 10
u_fourier = np.zeros_like(x)
for n in range(1, N_modes+1):
    # For initial condition u(x,0)=sin(πx), we have b1 = 1, others=0 ideally.
    # Here, for demonstration, include more modes with diminishing weight.
    b_n = 1/n if n==1 else 0.1/n  
    u_fourier += b_n * np.sin(n * np.pi * x / L) * np.exp(- (n * np.pi / L)**2 * a * t_heat)

### 2. Green's Function (Fundamental Solution) for Heat Equation
# ---------------------------------------------------------------
# The Green's function (heat kernel) for the 1D heat equation is given by:
#   G(x, t) = 1/(2√(π a t)) * exp(-x²/(4 a t))
# Here, we plot the response due to a point source centered at x0 = L/2.
def heat_green(x, t, a):
    return 1/(2*np.sqrt(np.pi * a * t)) * np.exp(-x**2/(4*a*t))

# Shift variable: compute G(x-L/2, t)
u_green = heat_green(x - L/2, t_heat, a)

### 3. Fourier Transform Method: Demonstration via the Heat Kernel
# ---------------------------------------------------------------
# For problems on infinite domain, Fourier transform of the heat kernel can be computed.
# Here, we compare numerical FFT of a discretized heat kernel against its analytical behavior.
# We use a spatial grid for x over a finite interval, and then compute FFT.
N_fft = 1024
L_fft = 10.0
dx_fft = L_fft / N_fft
x_fft = np.linspace(-L_fft/2, L_fft/2, N_fft, endpoint=False)
t_ft = t_heat  # use same t_heat as above
G_fft = heat_green(x_fft, t_ft, a)
# Compute FFT (shifted)
F_fft = fftshift(fft(G_fft))
freqs = fftshift(fftfreq(N_fft, d=dx_fft))
# Normalize for plotting amplitude
F_fft_norm = np.abs(F_fft) / N_fft

### 4. Special Function Expansion: Bessel Function J1(x)
# ---------------------------------------------------------------
# Bessel functions appear in PDEs in cylindrical coordinates.
# Here, we plot J1(x) as an example of special function expansion.
x_bessel = np.linspace(0, 20, 500)
J1 = jn(1, x_bessel)

### プロット（各セクションをまとめて表示） / Plot all sections
fig, axs = plt.subplots(2, 2, figsize=(14, 10))

# (1) Fourier series solution for Heat Equation
axs[0, 0].plot(x, u_fourier, 'r-', label=f"t = {t_heat}")
axs[0, 0].set_title("1. Heat Equation (Fourier Series)")
axs[0, 0].set_xlabel("x")
axs[0, 0].set_ylabel("u(x,t)")
axs[0, 0].grid(True)
axs[0, 0].legend()

# (2) Green's function solution for Heat Equation
axs[0, 1].plot(x, u_green, 'b-', label=f"t = {t_heat}")
axs[0, 1].set_title("2. Heat Kernel (Green's Function)")
axs[0, 1].set_xlabel("x")
axs[0, 1].set_ylabel("G(x-L/2, t)")
axs[0, 1].grid(True)
axs[0, 1].legend()

# (3) Fourier Transform (FFT) of Heat Kernel
axs[1, 0].plot(freqs, F_fft_norm, 'm-')
axs[1, 0].set_title("3. FFT of Heat Kernel")
axs[1, 0].set_xlabel("Frequency")
axs[1, 0].set_ylabel("|F(α)|")
axs[1, 0].grid(True)

# (4) Bessel Function J1(x) plot
axs[1, 1].plot(x_bessel, J1, 'g-')
axs[1, 1].set_title("4. Bessel Function J₁(x)")
axs[1, 1].set_xlabel("x")
axs[1, 1].set_ylabel("J₁(x)")
axs[1, 1].grid(True)

plt.tight_layout()
plt.show()

# プログラム名: heat_wave_fourier_solution_plots.py
# 内容: 添付画像に示された熱伝導方程式・波動方程式の解析解をPythonで可視化
# Purpose: Plot the exact solutions of heat and wave PDEs as shown in provided textbook images

import numpy as np
import matplotlib.pyplot as plt

# === 1次元熱伝導方程式の基本解（ガウス型）===
def u_heat_gauss(x, t):
    return (1 / np.sqrt(4 * np.pi * t)) * np.exp(-x**2 / (4 * t))

x_vals = np.linspace(-5, 5, 500)
t_vals = [0.1, 0.2, 0.4, 0.6, 0.8, 1.0]

plt.figure(figsize=(10, 6))
for t in t_vals:
    plt.plot(x_vals, u_heat_gauss(x_vals, t), label=f"t={t}")
plt.title("1D Heat Equation Gaussian Solution")
plt.xlabel("x")
plt.ylabel("u(x, t)")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# === フーリエ級数解 (熱伝導) from 画像2: u(x,t) = 20/pi * Σ 1/k * (1 - cos(kπ/2)) sin(kπx) e^{-k^2π^2 t} ===
def u_fourier_heat(x, t, K=60):
    u = 0
    for k in range(1, K+1):
        coeff = (1 - np.cos(k * np.pi / 2)) / k
        u += coeff * np.sin(k * np.pi * x) * np.exp(- (k * np.pi)**2 * t)
    return (20 / np.pi) * u

x_grid = np.linspace(0, 1, 500)
t_sample = 0.05
u_vals = u_fourier_heat(x_grid, t_sample)

plt.figure(figsize=(8, 4))
plt.plot(x_grid, u_vals)
plt.title(f"Fourier Series Solution of Heat Equation at t={t_sample}")
plt.xlabel("x")
plt.ylabel("u(x, t)")
plt.grid(True)
plt.tight_layout()
plt.show()

# === 2次元波動方程式のフーリエ解（画像3より）===
def u_2d_wave(x, y, t, N=10):
    result = 0
    for k in range(1, N+1):
        for j in range(1, N+1):
            lambda_kj = np.sqrt(k**2 + j**2) * np.pi / 4
            term = (1 / (k * j)) * np.sin(k * np.pi / 2) * np.sin(j * np.pi / 2)
            term *= np.sin(k * np.pi * x / 4) * np.sin(j * np.pi * y / 4)
            term *= np.cos(lambda_kj * t)
            result += term
    return (8 / np.pi**2) * result

# 2次元プロット（t=0）
x_2d = np.linspace(0, 4, 100)
y_2d = np.linspace(0, 4, 100)
X, Y = np.meshgrid(x_2d, y_2d)
Z = u_2d_wave(X, Y, t=0)

fig = plt.figure(figsize=(8, 5))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, cmap='viridis')
ax.set_title("2D Wave Equation Solution (t = 0)")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("u(x, y, 0)")
plt.tight_layout()
plt.show()

# プログラム名: pde_solutions_combined.py
# 内容: 偏微分方程式（熱伝導・波動・ラプラス方程式）の解析解または数値解を合体可視化
# Program: Combined visualization of heat, wave, and Laplace PDE solutions using analytical and numerical methods

import numpy as np
import matplotlib.pyplot as plt

# --- 1. 一次元熱伝導方程式（解析解）/ 1D Heat Equation (Analytical Solution) ---
def u_heat_1d(x, t):
    return (1 / np.sqrt(4 * np.pi * t)) * np.exp(-x**2 / (4 * t))

# --- 2. 一次元波動方程式（初期条件付き）/ 1D Wave Equation (Simple Analytical) ---
def u_wave_1d(x, t, L=1.0, c=1.0):
    return np.sin(np.pi * x / L) * np.cos(np.pi * c * t / L)

# --- 3. 二次元ラプラス方程式（簡易解析解）/ 2D Laplace Equation Solution ---
def u_laplace_2d(x, y, N=10):
    u = 0
    for n in range(1, N+1):
        coef = 4 / (n * np.pi) * (1 - (-1)**n)
        u += coef * np.sinh(n * np.pi * y) / np.sinh(n * np.pi) * np.sin(n * np.pi * x)
    return u

# --- グリッド定義 / Grid definition ---
x_1d = np.linspace(-5, 5, 500)
x_wave = np.linspace(0, 1, 500)
x2d = np.linspace(0, 1, 100)
y2d = np.linspace(0, 1, 100)
X2d, Y2d = np.meshgrid(x2d, y2d)

# --- プロット1: 一次元熱伝導方程式 ---
plt.figure(figsize=(12, 4))
for t in [0.05, 0.1, 0.2, 0.5, 1.0]:
    plt.plot(x_1d, u_heat_1d(x_1d, t), label=f"t={t}")
plt.title("1D Heat Equation (Analytical Solution)")
plt.xlabel("x")
plt.ylabel("u(x, t)")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# --- プロット2: 一次元波動方程式 ---
plt.figure(figsize=(12, 4))
for t in [0.0, 0.2, 0.4, 0.6, 0.8]:
    plt.plot(x_wave, u_wave_1d(x_wave, t), label=f"t={t}")
plt.title("1D Wave Equation (Analytical Solution)")
plt.xlabel("x")
plt.ylabel("u(x, t)")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# --- プロット3: 2D ラプラス方程式（定常解）---
Z = u_laplace_2d(X2d, Y2d)
fig = plt.figure(figsize=(8, 5))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X2d, Y2d, Z, cmap='viridis')
ax.set_title("2D Laplace Equation Solution")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("u(x, y)")
plt.tight_layout()
plt.show()

# プログラム名: pde_classification_examples.py
# 内容: 双曲型・放物型・楕円型の代表的なPDE（波動・熱・ラプラス）の可視化を統合
# Program Name: Combined plots for classification of PDEs (Hyperbolic, Parabolic, Elliptic)

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

# === 双曲型: 1D Wave Equation / Hyperbolic PDE ===
def u_wave_1d(x, t, c=1.0):
    return np.sin(np.pi * x) * np.cos(np.pi * c * t)

# === 放物型: 1D Heat Equation / Parabolic PDE ===
def u_heat_1d(x, t):
    return (1 / np.sqrt(4 * np.pi * t)) * np.exp(-x**2 / (4 * t))

# === 楕円型: 2D Laplace Equation / Elliptic PDE ===
def u_laplace_2d(x, y, N=10):
    u = 0
    for n in range(1, N+1):
        coeff = 4 / (n * np.pi) * (1 - (-1)**n)
        u += coeff * np.sinh(n * np.pi * y) / np.sinh(n * np.pi) * np.sin(n * np.pi * x)
    return u

# === 計算グリッド / Grids ===
x_1d = np.linspace(0, 1, 500)
x_2d = np.linspace(0, 1, 100)
y_2d = np.linspace(0, 1, 100)
X2d, Y2d = np.meshgrid(x_2d, y_2d)

# --- プロット1: 双曲型（波動） ---
plt.figure(figsize=(10, 4))
for t in [0.0, 0.2, 0.4, 0.6]:
    plt.plot(x_1d, u_wave_1d(x_1d, t), label=f"t={t}")
plt.title("Hyperbolic PDE (Wave Equation)")
plt.xlabel("x")
plt.ylabel("u(x,t)")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# --- プロット2: 放物型（熱） ---
plt.figure(figsize=(10, 4))
for t in [0.05, 0.1, 0.2, 0.5]:
    plt.plot(x_1d - 0.5, u_heat_1d(x_1d - 0.5, t), label=f"t={t}")
plt.title("Parabolic PDE (Heat Equation)")
plt.xlabel("x")
plt.ylabel("u(x,t)")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# --- プロット3: 楕円型（ラプラス） ---
Z = u_laplace_2d(X2d, Y2d)
fig = plt.figure(figsize=(8, 5))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X2d, Y2d, Z, cmap='plasma')
ax.set_title("Elliptic PDE (Laplace Equation)")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("u(x, y)")
plt.tight_layout()
plt.show()