PINN（Physics-Informed Neural Network、物理インフォームドニューラルネットワーク）は、ニューラルネットワークの学習に物理法則（偏微分方程式など）の制約を組み込む手法です。Raissi et al.（2019）によって体系化され、データ駆動アプローチと物理法則ベースアプローチの利点を融合する枠組みとして注目を集めています。
PINNは、観測データが少ない場合でも物理法則が正則化の役割を果たすため、汎化性能が向上します。また、偏微分方程式の数値解法（有限差分法、有限要素法など）の代替としても使え、メッシュ生成が不要という利点もあります。
本記事の内容
- PINNの動機と基本的な定式化
- 損失関数の構成（データ損失 + 物理損失 + 境界条件損失）
- 自動微分によるPDE残差の計算
- コロケーション点の配置戦略
- 境界条件・初期条件のハード/ソフト制約
- 逆問題への応用
- PyTorchによる1D熱方程式のPINN実装と解析解との比較
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
PINNの動機
偏微分方程式（PDE）の数値解法は科学技術計算の中核ですが、従来手法にはいくつかの課題があります。
- メッシュ依存性: 有限差分法や有限要素法は計算領域のメッシュ（格子）生成が必要で、複雑な形状では困難
- 高次元への拡張が困難: 3次元以上の問題では計算量が爆発的に増加
- 逆問題への適用が非自明: PDE係数の推定やデータ同化には追加の定式化が必要
PINNはこれらの課題に対して、ニューラルネットワークの万能近似性と自動微分の仕組みを利用して、統一的なフレームワークを提供します。
基本的な定式化
問題設定
領域 $\Omega \subset \mathbb{R}^d$ と時間区間 $[0, T]$ 上で定義された一般的なPDE問題を考えます。
$$ \begin{align} \mathcal{N}[u](\bm{x}, t) &= 0, \quad \bm{x} \in \Omega, \quad t \in [0, T] \\ \mathcal{B}[u](\bm{x}, t) &= g(\bm{x}, t), \quad \bm{x} \in \partial\Omega, \quad t \in [0, T] \\ u(\bm{x}, 0) &= u_0(\bm{x}), \quad \bm{x} \in \Omega \end{align} $$
ここで、
- $\mathcal{N}[\cdot]$: PDEの微分作用素
- $\mathcal{B}[\cdot]$: 境界条件の作用素
- $g$: 境界条件の関数
- $u_0$: 初期条件
ニューラルネットワーク近似
PINNでは、解 $u(\bm{x}, t)$ をパラメータ $\bm{\theta}$ を持つニューラルネットワーク $u_{\bm{\theta}}(\bm{x}, t)$ で近似します。
$$ u(\bm{x}, t) \approx u_{\bm{\theta}}(\bm{x}, t) = \text{NN}_{\bm{\theta}}(\bm{x}, t) $$
入力は空間座標 $\bm{x}$ と時刻 $t$、出力は解の近似値です。ネットワークは典型的にはフィードフォワードニューラルネットワーク（多層パーセプトロン）を使用します。
$$ u_{\bm{\theta}}(\bm{x}, t) = \bm{W}_L \sigma(\bm{W}_{L-1} \sigma(\cdots \sigma(\bm{W}_1 [\bm{x}; t] + \bm{b}_1) \cdots) + \bm{b}_{L-1}) + \bm{b}_L $$
ここで $\sigma$ は活性化関数（$\tanh$, $\sin$ など）、$[\bm{x}; t]$ は $\bm{x}$ と $t$ の結合です。
損失関数
PINNの損失関数は3つの項から構成されます。
$$ \mathcal{L}(\bm{\theta}) = \lambda_{\text{data}} \mathcal{L}_{\text{data}} + \lambda_{\text{physics}} \mathcal{L}_{\text{physics}} + \lambda_{\text{bc}} \mathcal{L}_{\text{bc}} $$
ここで $\lambda_{\text{data}}, \lambda_{\text{physics}}, \lambda_{\text{bc}}$ は各項の重み係数です。
データ損失: 観測データ $\{(\bm{x}_i^{\text{data}}, t_i^{\text{data}}, u_i^{\text{data}})\}_{i=1}^{N_{\text{data}}}$ に対する二乗誤差です。
$$ \mathcal{L}_{\text{data}} = \frac{1}{N_{\text{data}}} \sum_{i=1}^{N_{\text{data}}} \left| u_{\bm{\theta}}(\bm{x}_i^{\text{data}}, t_i^{\text{data}}) – u_i^{\text{data}} \right|^2 $$
物理損失（PDE残差損失）: 領域内のコロケーション点 $\{(\bm{x}_j^{\text{col}}, t_j^{\text{col}})\}_{j=1}^{N_{\text{col}}}$ におけるPDE残差の二乗誤差です。
$$ \mathcal{L}_{\text{physics}} = \frac{1}{N_{\text{col}}} \sum_{j=1}^{N_{\text{col}}} \left| \mathcal{N}[u_{\bm{\theta}}](\bm{x}_j^{\text{col}}, t_j^{\text{col}}) \right|^2 $$
PDE残差 $f(\bm{x}, t; \bm{\theta}) = \mathcal{N}[u_{\bm{\theta}}](\bm{x}, t)$ を計算するために、$u_{\bm{\theta}}$ の偏微分が必要です。これは自動微分によって正確に計算できます。
境界条件・初期条件損失: 境界と初期条件での誤差です。
$$ \mathcal{L}_{\text{bc}} = \frac{1}{N_{\text{bc}}} \sum_{k=1}^{N_{\text{bc}}} \left| \mathcal{B}[u_{\bm{\theta}}](\bm{x}_k^{\text{bc}}, t_k^{\text{bc}}) – g(\bm{x}_k^{\text{bc}}, t_k^{\text{bc}}) \right|^2 + \frac{1}{N_{\text{ic}}} \sum_{l=1}^{N_{\text{ic}}} \left| u_{\bm{\theta}}(\bm{x}_l^{\text{ic}}, 0) – u_0(\bm{x}_l^{\text{ic}}) \right|^2 $$
自動微分によるPDE残差の計算
PINNの核心的な仕組みは、ニューラルネットワークの出力に対する偏微分を自動微分（Automatic Differentiation）で正確に計算することです。
例えば、1次元の熱方程式 $\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}$ のPDE残差は以下です。
$$ f(x, t; \bm{\theta}) = \frac{\partial u_{\bm{\theta}}}{\partial t} – \alpha \frac{\partial^2 u_{\bm{\theta}}}{\partial x^2} $$
自動微分では、ネットワークの各層を通じた計算グラフをたどることで、以下の偏微分を正確に計算します。
$$ \frac{\partial u_{\bm{\theta}}}{\partial t} = \frac{\partial \text{NN}_{\bm{\theta}}}{\partial t}(x, t) $$
$$ \frac{\partial^2 u_{\bm{\theta}}}{\partial x^2} = \frac{\partial}{\partial x}\left(\frac{\partial u_{\bm{\theta}}}{\partial x}\right) $$
PyTorchでは torch.autograd.grad 関数を使って、入力変数に対する微分を計算できます。2階微分は、1階微分の結果に対してもう一度 torch.autograd.grad を適用することで得られます。
有限差分法による数値微分とは異なり、自動微分は浮動小数点精度の範囲で「正確」な微分を計算します。離散化誤差がないことがPINNの重要な特長です。
コロケーション点の配置戦略
物理損失を評価するコロケーション点 $\{(\bm{x}_j, t_j)\}$ の配置は、PINNの精度に大きく影響します。
一様ランダムサンプリング
最もシンプルな方法は、領域 $\Omega \times [0, T]$ 内で一様にランダムサンプリングすることです。
$$ \bm{x}_j \sim \text{Uniform}(\Omega), \quad t_j \sim \text{Uniform}([0, T]) $$
実装が容易で、多くの場合に良好に機能します。
Latin Hypercube Sampling（LHS）
LHSは各次元を $N$ 個の等間隔のストラタ（層）に分割し、各ストラタから1点ずつサンプリングする方法です。一様サンプリングよりも空間を均一にカバーでき、低次元問題では効果的です。
適応的サンプリング
学習の途中でPDE残差が大きい領域にコロケーション点を追加する戦略です。残差 $|f(\bm{x}, t; \bm{\theta})|$ が大きい領域は近似精度が低いため、追加点を配置することで効率的に精度を改善できます。
$$ p(\bm{x}, t) \propto |f(\bm{x}, t; \bm{\theta})|^k $$
この確率密度に従ってサンプリングすることで、難しい領域に計算リソースを集中させます。$k$ は集中度を制御するパラメータです。
境界条件のハード制約とソフト制約
ソフト制約
上述の損失関数に境界条件の項 $\mathcal{L}_{\text{bc}}$ を含める方法です。最も一般的なアプローチですが、境界条件が厳密には満たされず、重み $\lambda_{\text{bc}}$ の調整が必要です。
ハード制約
ネットワークの出力を変形して、境界条件を構造的に満たすようにする方法です。例えば、ディリクレ境界条件 $u(0, t) = 0$, $u(L, t) = 0$ の場合、以下のように構成します。
$$ u_{\bm{\theta}}(x, t) = x(L – x) \cdot \text{NN}_{\bm{\theta}}(x, t) $$
$x = 0$ および $x = L$ で自動的に $u = 0$ が満たされます。ネットワークが何を出力しても境界条件が厳密に満たされるため、学習が容易になります。
初期条件 $u(x, 0) = u_0(x)$ のハード制約は以下のように構成できます。
$$ u_{\bm{\theta}}(x, t) = u_0(x) + t \cdot \text{NN}_{\bm{\theta}}(x, t) $$
$t = 0$ で $u_{\bm{\theta}} = u_0(x)$ が保証されます。
ディリクレ境界条件と初期条件の両方を同時にハード制約する場合は、
$$ u_{\bm{\theta}}(x, t) = (1-t/T) \cdot u_0(x) + x(L-x) \cdot t \cdot \text{NN}_{\bm{\theta}}(x, t) $$
のように組み合わせます。ただし、複雑な境界条件ではハード制約の設計が困難な場合もあります。
逆問題への応用
PINNの強力な応用の一つが逆問題です。PDEの形は既知だが、一部の係数が未知の場合に、観測データからその係数を推定します。
例えば、熱方程式の熱拡散率 $\alpha$ が未知の場合を考えます。
$$ \frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2} $$
$\alpha$ を学習可能なパラメータとして、ネットワークのパラメータ $\bm{\theta}$ と同時に最適化します。
$$ \mathcal{L}(\bm{\theta}, \alpha) = \mathcal{L}_{\text{data}} + \lambda \mathcal{L}_{\text{physics}}(\bm{\theta}, \alpha) $$
ここでのポイントは、データ損失がPDE係数の推定に情報を提供し、物理損失がデータ間の補間に物理的な整合性を与えることです。この二つの損失の協調により、少数の観測データからでも $\alpha$ を正確に推定できる場合があります。
具体的には、$\alpha$ をPyTorchの nn.Parameter として宣言し、勾配降下法で $\bm{\theta}$ と $\alpha$ を同時に最適化します。
収束の難しさと解決策
PINNの学習は通常のニューラルネットワークに比べて収束が難しいことが知られています。
問題1: 損失項間の競合
データ損失と物理損失の勾配方向が異なる場合、一方を下げると他方が上がるという競合が生じます。重み $\lambda$ の調整が不適切だと、どちらか一方のみが最適化されます。
解決策: 学習の進行に応じて重みを動的に調整する方法が提案されています。例えば、各損失項の勾配のノルムを均等化するGradNorm法や、以下の適応的重み付けが使われます。
$$ \lambda_{\text{physics}}^{(k)} = \frac{\max_{\bm{\theta}} \|\nabla_{\bm{\theta}} \mathcal{L}_{\text{data}}\|}{\overline{\|\nabla_{\bm{\theta}} \mathcal{L}_{\text{physics}}\|}} $$
問題2: スペクトルバイアス
ニューラルネットワークは低周波成分を先に学習し、高周波成分の学習が遅いというスペクトルバイアスを持ちます。急激に変化するPDEの解（衝撃波など）では、この問題が顕著になります。
解決策: フーリエ特徴量を入力に加える方法や、活性化関数として $\sin$ を使うSIREN（Sinusoidal Representation Networks）が効果的です。
$$ \bm{\gamma}(\bm{x}) = [\sin(2\pi \bm{B} \bm{x}), \cos(2\pi \bm{B} \bm{x})] $$
ここで $\bm{B}$ はランダムに初期化されたフーリエ基底行列です。
問題3: 学習率と最適化
PINNの損失関数は高度に非凸であり、局所最適解に陥りやすいです。
解決策:
- Adamオプティマイザで初期の大域的探索を行い、L-BFGSで微調整する2段階最適化
- 学習率のウォームアップとコサインアニーリング
- 残差接続やレイヤー正規化によるネットワーク設計の工夫
Pythonでの実装: 1D熱方程式
ここからは、1次元熱方程式のPINNをPyTorchで実装します。
問題設定
以下の1次元熱方程式を考えます。
$$ \frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}, \quad x \in [0, L], \quad t \in [0, T] $$
境界条件: $u(0, t) = 0$, $u(L, t) = 0$
初期条件: $u(x, 0) = \sin(\pi x / L)$
この問題の解析解は以下です。
$$ u(x, t) = \sin\left(\frac{\pi x}{L}\right) \exp\left(-\alpha \left(\frac{\pi}{L}\right)^2 t\right) $$
導出を確認しましょう。変数分離法で $u(x, t) = X(x) T(t)$ と置くと、熱方程式は以下になります。
$$ X(x) T'(t) = \alpha X”(x) T(t) $$
両辺を $\alpha X(x) T(t)$ で割ると、
$$ \frac{T'(t)}{\alpha T(t)} = \frac{X”(x)}{X(x)} = -\lambda $$
ここで $\lambda$ は分離定数です。右辺の空間部分は境界条件 $X(0) = X(L) = 0$ のもとで解くと、
$$ X”(x) + \lambda X(x) = 0 \implies X_n(x) = \sin\left(\frac{n\pi x}{L}\right), \quad \lambda_n = \left(\frac{n\pi}{L}\right)^2 $$
時間部分は、
$$ T'(t) = -\alpha \lambda_n T(t) \implies T_n(t) = \exp\left(-\alpha \left(\frac{n\pi}{L}\right)^2 t\right) $$
初期条件 $u(x, 0) = \sin(\pi x / L)$ は $n=1$ のモードのみを含むので、解析解は上記の通りです。
実装
import torch
import torch.nn as nn
import numpy as np
import matplotlib.pyplot as plt
# 再現性
torch.manual_seed(42)
np.random.seed(42)
# 物理パラメータ
alpha_true = 0.01 # 熱拡散率
L = 1.0 # 空間領域の長さ
T_final = 1.0 # 最終時刻
# 解析解
def analytical_solution(x, t, alpha=alpha_true, L_val=L):
return np.sin(np.pi * x / L_val) * np.exp(-alpha * (np.pi / L_val)**2 * t)
# デバイス設定
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
ニューラルネットワークの定義
class PINN(nn.Module):
"""Physics-Informed Neural Network for 1D Heat Equation"""
def __init__(self, hidden_dim=64, num_layers=4):
super().__init__()
layers = []
# 入力層: (x, t) -> hidden_dim
layers.append(nn.Linear(2, hidden_dim))
layers.append(nn.Tanh())
# 隠れ層
for _ in range(num_layers - 1):
layers.append(nn.Linear(hidden_dim, hidden_dim))
layers.append(nn.Tanh())
# 出力層: hidden_dim -> 1
layers.append(nn.Linear(hidden_dim, 1))
self.network = nn.Sequential(*layers)
# 重みの初期化（Xavier初期化）
for m in self.network:
if isinstance(m, nn.Linear):
nn.init.xavier_normal_(m.weight)
nn.init.zeros_(m.bias)
def forward(self, x, t):
"""順伝播: (x, t) -> u"""
inputs = torch.cat([x, t], dim=-1)
return self.network(inputs)
PDE残差の計算
def compute_pde_residual(model, x, t, alpha=alpha_true):
"""
自動微分によるPDE残差の計算
f = u_t - alpha * u_xx = 0
"""
# 勾配計算を有効にする
x.requires_grad_(True)
t.requires_grad_(True)
u = model(x, t)
# 1階微分の計算
u_t = torch.autograd.grad(
u, t,
grad_outputs=torch.ones_like(u),
create_graph=True, retain_graph=True
)[0]
u_x = torch.autograd.grad(
u, x,
grad_outputs=torch.ones_like(u),
create_graph=True, retain_graph=True
)[0]
# 2階微分の計算
u_xx = torch.autograd.grad(
u_x, x,
grad_outputs=torch.ones_like(u_x),
create_graph=True, retain_graph=True
)[0]
# PDE残差: u_t - alpha * u_xx
residual = u_t - alpha * u_xx
return residual
コロケーション点の生成と学習
def generate_training_points(n_collocation=5000, n_bc=200, n_ic=200):
"""訓練用のコロケーション点・境界条件点・初期条件点を生成"""
# コロケーション点（領域内部）
x_col = torch.rand(n_collocation, 1) * L
t_col = torch.rand(n_collocation, 1) * T_final
# 境界条件点: x=0 と x=L
t_bc = torch.rand(n_bc, 1) * T_final
x_bc_left = torch.zeros(n_bc, 1) # x = 0
x_bc_right = torch.ones(n_bc, 1) * L # x = L
u_bc = torch.zeros(n_bc, 1) # u = 0 at boundaries
# 初期条件点: t=0
x_ic = torch.rand(n_ic, 1) * L
t_ic = torch.zeros(n_ic, 1)
u_ic = torch.sin(np.pi * x_ic / L) # u(x,0) = sin(πx/L)
return (x_col, t_col,
x_bc_left, x_bc_right, t_bc, u_bc,
x_ic, t_ic, u_ic)
# モデルの初期化
model = PINN(hidden_dim=64, num_layers=4).to(device)
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
scheduler = torch.optim.lr_scheduler.StepLR(optimizer, step_size=2000, gamma=0.5)
# 訓練データの生成
(x_col, t_col,
x_bc_left, x_bc_right, t_bc, u_bc,
x_ic, t_ic, u_ic) = generate_training_points()
# デバイスに転送
x_col, t_col = x_col.to(device), t_col.to(device)
x_bc_left, x_bc_right = x_bc_left.to(device), x_bc_right.to(device)
t_bc, u_bc = t_bc.to(device), u_bc.to(device)
x_ic, t_ic, u_ic = x_ic.to(device), t_ic.to(device), u_ic.to(device)
# 損失の重み
lambda_physics = 1.0
lambda_bc = 10.0
lambda_ic = 10.0
# 学習ループ
epochs = 10000
loss_history = {'total': [], 'physics': [], 'bc': [], 'ic': []}
for epoch in range(1, epochs + 1):
model.train()
optimizer.zero_grad()
# 物理損失（PDE残差）
residual = compute_pde_residual(model, x_col, t_col)
loss_physics = torch.mean(residual ** 2)
# 境界条件損失
u_pred_left = model(x_bc_left, t_bc)
u_pred_right = model(x_bc_right, t_bc)
loss_bc = torch.mean((u_pred_left - u_bc)**2) + torch.mean((u_pred_right - u_bc)**2)
# 初期条件損失
u_pred_ic = model(x_ic, t_ic)
loss_ic = torch.mean((u_pred_ic - u_ic)**2)
# 総損失
loss = lambda_physics * loss_physics + lambda_bc * loss_bc + lambda_ic * loss_ic
loss.backward()
optimizer.step()
scheduler.step()
# 記録
loss_history['total'].append(loss.item())
loss_history['physics'].append(loss_physics.item())
loss_history['bc'].append(loss_bc.item())
loss_history['ic'].append(loss_ic.item())
if epoch % 2000 == 0:
print(f'Epoch {epoch}/{epochs}, '
f'Total: {loss.item():.6f}, '
f'Physics: {loss_physics.item():.6f}, '
f'BC: {loss_bc.item():.6f}, '
f'IC: {loss_ic.item():.6f}')
学習曲線の可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 総損失
axes[0].semilogy(loss_history['total'], linewidth=1)
axes[0].set_xlabel('Epoch', fontsize=12)
axes[0].set_ylabel('Loss', fontsize=12)
axes[0].set_title('Total Loss', fontsize=14)
axes[0].grid(True, alpha=0.3)
# 各項の損失
axes[1].semilogy(loss_history['physics'], label='Physics', linewidth=1)
axes[1].semilogy(loss_history['bc'], label='BC', linewidth=1)
axes[1].semilogy(loss_history['ic'], label='IC', linewidth=1)
axes[1].set_xlabel('Epoch', fontsize=12)
axes[1].set_ylabel('Loss', fontsize=12)
axes[1].set_title('Loss Components', fontsize=14)
axes[1].legend(fontsize=10)
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
解析解との比較
# 評価用の格子点
nx, nt = 100, 100
x_eval = np.linspace(0, L, nx)
t_eval = np.linspace(0, T_final, nt)
X_mesh, T_mesh = np.meshgrid(x_eval, t_eval)
# PINN予測
model.eval()
x_flat = torch.tensor(X_mesh.flatten(), dtype=torch.float32).unsqueeze(1).to(device)
t_flat = torch.tensor(T_mesh.flatten(), dtype=torch.float32).unsqueeze(1).to(device)
with torch.no_grad():
u_pred = model(x_flat, t_flat).cpu().numpy().reshape(nt, nx)
# 解析解
u_exact = analytical_solution(X_mesh, T_mesh)
# 誤差
error = np.abs(u_pred - u_exact)
print(f'Mean Absolute Error: {np.mean(error):.6f}')
print(f'Max Absolute Error: {np.max(error):.6f}')
print(f'Relative L2 Error: {np.linalg.norm(u_pred - u_exact) / np.linalg.norm(u_exact):.6f}')
# 可視化
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
# PINN予測
im0 = axes[0].contourf(X_mesh, T_mesh, u_pred, levels=50, cmap='viridis')
axes[0].set_xlabel('x', fontsize=12)
axes[0].set_ylabel('t', fontsize=12)
axes[0].set_title('PINN Prediction', fontsize=14)
plt.colorbar(im0, ax=axes[0])
# 解析解
im1 = axes[1].contourf(X_mesh, T_mesh, u_exact, levels=50, cmap='viridis')
axes[1].set_xlabel('x', fontsize=12)
axes[1].set_ylabel('t', fontsize=12)
axes[1].set_title('Analytical Solution', fontsize=14)
plt.colorbar(im1, ax=axes[1])
# 絶対誤差
im2 = axes[2].contourf(X_mesh, T_mesh, error, levels=50, cmap='hot')
axes[2].set_xlabel('x', fontsize=12)
axes[2].set_ylabel('t', fontsize=12)
axes[2].set_title('Absolute Error', fontsize=14)
plt.colorbar(im2, ax=axes[2])
plt.tight_layout()
plt.show()
特定時刻での断面比較
# 特定時刻でのu(x)の比較
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
time_snapshots = [0.0, 0.3, 0.7]
for ax, t_snap in zip(axes, time_snapshots):
# PINN予測
x_plot = torch.tensor(x_eval, dtype=torch.float32).unsqueeze(1).to(device)
t_plot = torch.full_like(x_plot, t_snap)
with torch.no_grad():
u_pinn = model(x_plot, t_plot).cpu().numpy().flatten()
# 解析解
u_anal = analytical_solution(x_eval, t_snap)
ax.plot(x_eval, u_anal, 'b-', linewidth=2, label='Analytical')
ax.plot(x_eval, u_pinn, 'r--', linewidth=2, label='PINN')
ax.set_xlabel('x', fontsize=12)
ax.set_ylabel('u(x, t)', fontsize=12)
ax.set_title(f't = {t_snap:.1f}', fontsize=14)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
plt.suptitle('PINN vs Analytical Solution at Different Times', fontsize=14)
plt.tight_layout()
plt.show()
逆問題: 熱拡散率の推定
class PINN_Inverse(nn.Module):
"""逆問題用PINN: 熱拡散率αを同時に推定"""
def __init__(self, hidden_dim=64, num_layers=4):
super().__init__()
layers = []
layers.append(nn.Linear(2, hidden_dim))
layers.append(nn.Tanh())
for _ in range(num_layers - 1):
layers.append(nn.Linear(hidden_dim, hidden_dim))
layers.append(nn.Tanh())
layers.append(nn.Linear(hidden_dim, 1))
self.network = nn.Sequential(*layers)
for m in self.network:
if isinstance(m, nn.Linear):
nn.init.xavier_normal_(m.weight)
nn.init.zeros_(m.bias)
# 推定するパラメータ: 熱拡散率α（初期値は真の値から離れた値）
self.log_alpha = nn.Parameter(torch.tensor(np.log(0.1))) # 初期値0.1（真の値は0.01）
@property
def alpha(self):
return torch.exp(self.log_alpha) # 正の値を保証
def forward(self, x, t):
inputs = torch.cat([x, t], dim=-1)
return self.network(inputs)
# 観測データの生成（解析解からサンプリング + ノイズ）
n_obs = 200
x_obs = np.random.rand(n_obs) * L
t_obs = np.random.rand(n_obs) * T_final
u_obs = analytical_solution(x_obs, t_obs) + 0.01 * np.random.randn(n_obs)
x_obs_t = torch.tensor(x_obs, dtype=torch.float32).unsqueeze(1).to(device)
t_obs_t = torch.tensor(t_obs, dtype=torch.float32).unsqueeze(1).to(device)
u_obs_t = torch.tensor(u_obs, dtype=torch.float32).unsqueeze(1).to(device)
# 逆問題モデルの学習
model_inv = PINN_Inverse(hidden_dim=64, num_layers=4).to(device)
optimizer_inv = torch.optim.Adam(model_inv.parameters(), lr=1e-3)
alpha_history = []
epochs_inv = 8000
for epoch in range(1, epochs_inv + 1):
model_inv.train()
optimizer_inv.zero_grad()
# データ損失
u_pred_data = model_inv(x_obs_t, t_obs_t)
loss_data = torch.mean((u_pred_data - u_obs_t)**2)
# PDE残差損失（αはモデルから取得）
x_col_inv = (torch.rand(3000, 1) * L).to(device)
t_col_inv = (torch.rand(3000, 1) * T_final).to(device)
x_col_inv.requires_grad_(True)
t_col_inv.requires_grad_(True)
u_col = model_inv(x_col_inv, t_col_inv)
u_t = torch.autograd.grad(u_col, t_col_inv, torch.ones_like(u_col),
create_graph=True)[0]
u_x = torch.autograd.grad(u_col, x_col_inv, torch.ones_like(u_col),
create_graph=True)[0]
u_xx = torch.autograd.grad(u_x, x_col_inv, torch.ones_like(u_x),
create_graph=True)[0]
residual_inv = u_t - model_inv.alpha * u_xx
loss_physics_inv = torch.mean(residual_inv**2)
loss_inv = loss_data + loss_physics_inv
loss_inv.backward()
optimizer_inv.step()
alpha_history.append(model_inv.alpha.item())
if epoch % 2000 == 0:
print(f'Epoch {epoch}/{epochs_inv}, '
f'Loss: {loss_inv.item():.6f}, '
f'Estimated α: {model_inv.alpha.item():.6f} '
f'(True: {alpha_true})')
# αの推定過程の可視化
plt.figure(figsize=(8, 5))
plt.plot(alpha_history, 'b-', linewidth=1.5, label='Estimated $\\alpha$')
plt.axhline(y=alpha_true, color='r', linestyle='--', linewidth=2, label=f'True $\\alpha$ = {alpha_true}')
plt.xlabel('Epoch', fontsize=12)
plt.ylabel('$\\alpha$', fontsize=14)
plt.title('Inverse Problem: Estimation of Thermal Diffusivity', fontsize=14)
plt.legend(fontsize=12)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print(f'Final estimated α: {alpha_history[-1]:.6f}')
print(f'True α: {alpha_true}')
print(f'Relative error: {abs(alpha_history[-1] - alpha_true) / alpha_true * 100:.2f}%')
まとめ
本記事では、PINN（Physics-Informed Neural Network）の理論と実装を解説しました。
- PINNはニューラルネットワークの損失関数にPDE残差 $\|\mathcal{N}[u_{\bm{\theta}}]\|^2$ を組み込むことで、物理法則を学習に反映する
- 自動微分により、PDE残差の計算が正確かつ効率的に行える。有限差分法と異なり離散化誤差がない
- コロケーション点の配置は一様ランダムサンプリングが基本だが、適応的サンプリングで効率化できる
- 境界条件はソフト制約（損失関数に追加）とハード制約（出力を構造的に変形）の2つのアプローチがある
- 逆問題（PDE係数の推定）にも自然に拡張でき、観測データと物理法則の協調により少数データでの推定が可能
- 1D熱方程式のPINN実装では、解析解との良好な一致と、逆問題での熱拡散率の推定を確認した
次のステップとして、以下の記事も参考にしてください。