微分方程式は、物理現象のモデル化から金融工学まで、科学技術の幅広い分野で不可欠なツールです。しかし、その解析は容易ではなく、多くの場合、複雑な数値計算が必要とされます。近年、この古典的な問題に対して、ニューラルネットワークを用いた新しいアプローチが注目を集めています。

今回は、PyTorchベースのライブラリである neurodiffeq を使い、ニューラルネットワークで微分方程式を解く方法を紹介します。neurodiffeq は、微分方程式の解をニューラルネットワークで近似するための強力かつ柔軟なフレームワークです。

この記事を通して、neurodiffeq の基本的なインストール方法から、簡単な常微分方程式系（ODEシステム）を例題として、実際にコードの実装例を紹介します。ニューラルネットワークがどのように微分方程式の解を学習していくのか、その一端に触れてみましょう。

微分方程式とは？

微分方程式

微分方程式とは、一言で言えば「未知の関数とその導関数（変化率）の関係を表す方程式」です。私たちの身の回りにある現象の多くは、時間や空間における変化として捉えることができ、微分方程式はこれらの変化の法則を数学的に記述するための強力な言語となります。

例えば、物体の運動（ニュートンの運動方程式）、熱の伝わり方（熱伝導方程式）、化学反応の速度、人口の増減モデル、感染症の拡大モデル、さらには経済現象の予測など、多くの事象を微分方程式として表現できます。これらの現象を理解し、予測するためには、これらの微分方程式を解くことが不可欠です。

微分方程式は、大きく分けて常微分方程式 (Ordinary Differential Equations, ODEs) と偏微分方程式 (Partial Differential Equations, PDEs) の2種類があります。

• 常微分方程式 (ODEs): 1つの独立変数（例えば時間）に関する未知関数とその導関数を含む方程式です。例えば、単振動のモデル \(\displaystyle m \frac{d^2x}{dt^2}​=−kx\) はODEです。
• 偏微分方程式 (PDEs): 2つ以上の独立変数に関する未知関数とその偏導関数を含む方程式です。例えば、弦の振動を表す波動方程式 \(\displaystyle \frac{\partial^2 u}{ \partial t^2}​=c^2 \frac{\partial^ 2 u}{ \partial x^2}\)​ はPDEです。

従来の微分方程式の解法

微分方程式を解く伝統的なアプローチには、解析解法と数値解法があります。

• 解析解法: 数式を変形したり、特定の関数形を仮定したりすることで、厳密な解を求める方法です。しかし、解析的に解ける微分方程式は限られており、多くの実用的な問題では解析解を得ることは困難です。
• 数値解法: 解析解が得られない場合に、コンピュータを用いて近似的な解を求める方法です。代表的なものに、オイラー法、ルンゲ＝クッタ法（いずれもODE向け）、有限差分法、有限要素法、スペクトル法（主にPDE向け）などがあります。これらの手法は、定義域を離散的な点（格子点）に分割し、各点での関数の値を反復計算などによって近似的に求めます。

これらの従来の数値解法は非常に強力で、多くの問題で有効性が示されていますが、いくつかの課題も抱えています。例えば、高次元の問題になると計算コストが爆発的に増加する「次元の呪い」や、複雑な形状の境界を持つ問題への対応、解の関数形が特殊な場合に精度が低下する可能性などが挙げられます。

なぜ、ニューラルネットワークで微分方程式を解くのか？

近年、これらの課題に対する新しいアプローチとして、ニューラルネットワーク（特に深層学習）を用いて微分方程式を解く研究が活発に進められています。このアプローチには、いくつかの魅力的な利点があります。

1. 連続的で微分可能な解の表現: 従来の数値解法の多くは離散点での解を与えるのに対し、ニューラルネットワークは入力変数から出力（解の値）を計算する連続的な関数を学習します。学習されたネットワークは微分可能であるため、解の導関数も容易に得られます。これは、物理法則に基づく制約を損失関数に組み込む際などに有利に働きます。
2. 複雑な関数の近似能力: ニューラルネットワークは「普遍性定理 (Universal Approximation Theorem)」により、十分な数のニューロンと適切な活性化関数を用いれば、任意の連続関数を任意の精度で近似できることが知られています。この性質により、複雑な形状の解や、従来の基底関数では表現しにくい解も捉えられる可能性があります。
3. 高次元問題への潜在的なスケーラビリティ: ニューラルネットワークを用いた手法（特にPhysics-Informed Neural Networks, PINNsなど）は、サンプリングベースで損失を計算するため、原理的にはメッシュ生成が不要であり、高次元問題に対しても従来のメッシュベースの手法よりスケーラブルである可能性が期待されています。
4. 逆問題への応用: 微分方程式のパラメータを観測データから推定する「逆問題」に対しても、ニューラルネットワークは有望な枠組みを提供します。微分方程式の項と観測データの両方を損失関数に組み込むことで、パラメータ推定と解の発見を同時に行うことができます。

これらの理由から、ニューラルネットワークは微分方程式を解くための新しい強力なツールとして期待されており、neurodiffeq のようなライブラリは、このアプローチをより多くの研究者やエンジニアが利用できるようにすることを目指しています。

neurodiffeq とは何か？

neurodiffeq は、ニューラルネットワークを用いて微分方程式を解くために開発された、Pythonのオープンソースライブラリです。その最大の特徴は、深層学習フレームワークである PyTorch を基盤としている点にあります。これにより、PyTorchが持つ自動微分機能や豊富なニューラルネットワークモジュール、最適化アルゴリズムといった強力なエコシステムを最大限に活用できます。

このライブラリは、常微分方程式 (ODEs) だけでなく、偏微分方程式 (PDEs) も対象としており、単一の方程式から連立方程式まで、幅広い種類の問題に対応可能です。実際に、neurodiffeq は学術研究の現場でも利用されており、例えばハーバード大学の Institute for Applied Computational Science (IACS) など、世界中の複数の研究グループでその有効性が試されています。

従来の数値解法、例えば有限差分法や有限要素法が、定義域を離散化し、各格子点での近似値を求めるのに対し、neurodiffeq は異なるアプローチを取ります。それは、解そのものをニューラルネットワークという連続関数で表現するというものです。学習が完了すると、得られるのは単なる離散的な数値の集合ではなく、定義域内の任意の点で値を評価でき、かつ微分可能な関数形式の解です。これは、解の物理的な解釈や、さらなる解析を行う上で大きな利点となります。

このアプローチの根底には、人工ニューラルネットワーク (ANN) が持つ「普遍関数近似能力」があります。これは、十分な表現力を持つANNは、任意の連続関数を望む精度で近似できるという理論です。この性質を利用し、微分方程式の解という未知の関数をANNで表現し、微分方程式そのものと初期条件・境界条件を満たすようにネットワークのパラメータ（重みやバイアス）を学習させるのです。

neurodiffeq の開発目標は、このANNを用いた微分方程式解法に関する既存の技術やアイデアを、より多くのユーザーが自身の問題に合わせて柔軟に利用できる形で実装し、提供することにあります。単に特定の問題を解くだけでなく、新しいアイデアを試したり、複雑な条件を設定したりするための使いやすいインターフェースとツール群を提供することで、この分野の研究と応用の発展に貢献することを目指しています。

neurodiffeq を使ってみよう！

このセクションでは、実際に neurodiffeq を使って微分方程式を解く手順を解説します。まずはインストール方法から見ていきましょう。

PyTorch をインストールしておくことを推奨します。これにより、環境に合わせた PyTorch のバージョン（CPU版、CUDA版など）を確実に指定できます。

PyTorch のインストールは、PyTorch公式サイト の指示に従って行うのが最も確実です。例えば、GPU版（CUDA 12.6）の PyTorch をインストールする場合は、以下のようなコマンドになります（バージョンは適宜確認してください）。

$ pip install torch torchvision torchaudio

PyTorch がインストールされたら、neurodiffeq をインストールします。-U オプションを付けることで、既にインストールされている場合は最新版にアップグレードされます。

$ pip install -U neurodiffeq

これで準備は完了です。より詳細なインストールオプションや、開発版のインストール方法については、neurodiffeq の公式ドキュメント を参照してください。

基本的な使用例 (ODEシステム：Lotka–Volterra方程式)

ここでは、捕食者と被食者の個体数変動を表す古典的なモデルである Lotka–Volterra方程式 を例に、neurodiffeq の基本的な使い方を説明します。

Lotka–Volterra方程式は、以下のような連立常微分方程式で表されます。

• 被食者 (\(x\)) の増殖率: \(\displaystyle \frac{dx}{dt}​=\alpha x − \beta x y\)
• 捕食者 (\(y\)) の増殖率: \(\displaystyle \frac{dy}{dt}​ = \delta x y − \gamma y\)

ここで、\(t\) は時間、\(x(t)\) は被食者の個体数、\(y(t)\) は捕食者の個体数を表します。\(\alpha\),\(\beta\),\(\gamma\),\(\delta\) は正の定数で、それぞれ被食者の自然増殖率、捕食による被食者の減少率、捕食による捕食者の増加率、捕食者の自然死亡率を意味します。

以下に、neurodiffeq を用いてこの方程式を解く手順を示します。

まず、必要なモジュールをインポートします。

import torch
import numpy as np
from neurodiffeq.neurodiffeq import diff
from neurodiffeq.solvers import Solver1D
from neurodiffeq.conditions import IVP
from neurodiffeq.networks import FCNN, SinActv
import matplotlib.pyplot as plt

次に、解きたい微分方程式をPythonの関数として定義します。neurodiffeq では、diff(関数, 独立変数, 微分の階数) という形で導関数を表現します。なお、方程式が0になる形（残差形式）でリストとして返す必要があります。

# パラメータの設定 (例)
alpha = 1.0
beta = 1.0
delta = 1.0
gamma = 1.0

def ode_system(x, y, t):
    # x, y は t の関数として扱われる
    # diff(x, t) は dx/dt を意味する
    # diff(y, t) は dy/dt を意味する
    return [
        diff(x, t) - (alpha * x - beta * x * y),
        diff(y, t) - (delta * x * y - gamma * y)
    ]

次に初期条件を定義します。初期値問題 (Initial Value Problem, IVP) の場合、初期条件は必須です。ここでは、\(t=0\)における \(x\)と \(y\)の値を設定します。具体的には、IVP オブジェクトのリストとして定義し、各要素がそれぞれの方程式の解（この場合は \(x(t)\) と \(y(t)\)）に対応します。

# 初期条件 (例: t=0 で x=1.5, y=1.0)
conditions = [
    IVP(t_0=0.0, u_0=1.5), # x(0) = 1.5
    IVP(t_0=0.0, u_0=1.0)  # y(0) = 1.0
]

次に、解関数を近似するためのニューラルネットワークを定義します。neurodiffeq はデフォルトで FCNN (Fully Connected Neural Network) を提供しています。活性化関数には SinActv (サイン関数をベースにした活性化関数) を使うと、周期的な解の表現に適している場合があります。

# 各解関数 (x(t), y(t)) に対するネットワークをリストで指定
nets = [
    FCNN(actv=SinActv),  # x(t)用
    FCNN(actv=SinActv),  # y(t)用
]

ここまでできたら、Solver1D (1次元問題用ソルバー) を設定し、学習を実行します。

# 解を求めたい時間範囲
t_min = 0.1
t_max = 12.0 # 例えば12単位時間まで

# ソルバーのインスタンス化
solver = Solver1D(
    ode_system=ode_system,
    conditions=conditions,
    t_min=t_min,
    t_max=t_max,
    nets=nets
)

# 学習の実行 (max_epochs で学習エポック数を指定)
solver.fit(max_epochs=3000)

学習が完了したら、get_solution() メソッドで解を取得できます。このメソッドは、callable なオブジェクト（関数のように呼び出せるオブジェクト）を返します。

# 解の取得
solution = solver.get_solution()

# 特定の時間点での解を評価
t_eval = np.linspace(t_min, t_max, 100) # 評価する時間点
x_sol, y_sol = solution(t_eval, to_numpy=True) # to_numpy=True でNumPy配列として取得

# 結果の表示 (例)
for i in range(5):
    print(f"t = {t_eval[i*20]:.2f}: x(t) = {x_sol[i*20]:.4f}, y(t) = {y_sol[i*20]:.4f}")

t = 0.10: x(t) = 1.4965, y(t) = 1.0511
t = 2.50: x(t) = 0.6496, y(t) = 1.1542
t = 4.91: x(t) = 1.0433, y(t) = 0.6376
t = 7.31: x(t) = 1.1929, y(t) = 1.4223
t = 9.72: x(t) = 0.6601, y(t) = 0.8581

最後に、得られた解をプロットして、個体数の時間変化や位相空間での軌跡を確認してみましょう。

plt.figure(figsize=(12, 5))

plt.subplot(1, 2, 1)
plt.plot(t_eval, x_sol, label='x(t) - Preys')
plt.plot(t_eval, y_sol, label='y(t) - Predators')
plt.xlabel('Time (t)')
plt.ylabel('Population')
plt.title('Solution of the Lotka-Volterra equation (time evolution)')
plt.legend()
plt.grid(True)

plt.subplot(1, 2, 2)
plt.plot(x_sol, y_sol)
plt.xlabel('Prey population x(t)')
plt.ylabel('Predator population y(t)')
plt.title('Phase space trajectories (x-y plot)')
plt.grid(True)

plt.tight_layout()
plt.show()

このように、neurodiffeq を使うと、微分方程式の定義、条件設定、ネットワーク構築、そして学習と解の評価までを比較的直感的に行うことができます。

より高度な話題

neurodiffeq は基本的な使い方以外にも、より複雑な問題設定や詳細な学習制御に対応するための様々な機能を提供しています。ここでは、その主要なコンポーネントと高度な機能について紹介します。

neurodiffeq の主要コンポーネント

neurodiffeq を使った微分方程式の求解は、主に以下の3つのコンポーネントを組み合わせて行われます。

• ソルバー (Solvers): 微分方程式の種類や解きたい問題の次元に応じて適切なソルバーを選択します。
  • Solver1D: 1次元の常微分方程式 (ODE) や偏微分方程式 (PDE) 用の基本的なソルバーです。上の例で使用しました。
  • Solver2D: 2次元のPDE用ソルバーです。例えば、時間 t と空間 x の2つの独立変数を持つ問題に使用します。
  • BundleSolver1D: 複数のパラメータセットを持つ微分方程式群の解を一度に効率よく求めるためのソルバーです。 他にも、問題の特性に応じたソルバーが用意されています。
• 条件 (Conditions): 微分方程式の解を一意に定めるためには、初期条件や境界条件が必要です。neurodiffeq では、これらの条件をオブジェクトとして定義し、ソルバーに渡します。
  • IVP (Initial Value Problem): 常微分方程式の初期値問題の条件を指定します。上の例で使用しました。
  • DirichletBVP2D (Dirichlet Boundary Value Problem for 2D): 2次元問題におけるディリクレ境界条件（境界上での関数の値を指定）を定義します。
  • NeumannBVP2D: ノイマン境界条件（境界上での関数の法線微分の値を指定）を定義します。
  • BundleIVP: BundleSolver で使用する、パラメータ化された初期条件を定義します。 これらの他にも、周期境界条件など、様々な種類の条件クラスが提供されています。
• ネットワーク (Networks): 解関数を近似するためのニューラルネットワークを指定します。
  • FCNN (Fully Connected Neural Network): neurodiffeq にデフォルトで用意されている全結合ニューラルネットワークです。隠れ層の数やユニット数、活性化関数などを指定できます。
  • カスタム torch.nn.Module: PyTorchユーザーであれば、より複雑な、あるいは問題に特化した独自のネットワークアーキテクチャを torch.nn.Module のサブクラスとして定義し、それを neurodiffeq のソルバーに渡すことができます。これにより、畳み込みニューラルネットワーク (CNN) やリカレントニューラルネットワーク (RNN) など、様々なネットワーク構造を利用する道が開かれます。 重要な点として、ソルバーには 解くべき未知関数の数だけネットワークのリストを渡す 必要があります。例えば、Lotka-Volterra方程式では \(x(t)\) と \(y(t)\) の2つの未知関数があったため、2つのネットワークをリスト nets として渡しました。

学習の制御と可視化

ニューラルネットワークの学習プロセスを適切に制御し、結果を評価することは非常に重要です。neurodiffeq はそのためのツールも提供しています。

• Monitor (モニター): 学習中の損失関数の値や、解の挙動をリアルタイムで可視化するための強力なツールです。特に Jupyter Notebook や Jupyter Lab 環境では、インタラクティブなプロットを生成し、学習がうまく進んでいるか、あるいは発散していないかなどを視覚的に確認できます。
• Sampling Strategies (Generators – ジェネレータ): ニューラルネットワークの学習では、微分方程式の残差や境界条件の誤差を計算するために、定義域内から点をサンプリングします。このサンプリング方法（点の数、分布、範囲など）は、学習の効率や解の精度に影響を与えることがあります。neurodiffeq では、このサンプリング戦略を柔軟に制御するためのジェネレータ機能を提供しています。
• Generator1D, Generator2D: それぞれ1次元、2次元問題用の基本的なサンプリングジェネレータです。一様分布や特定の範囲に重点を置いた分布など、様々なサンプリング方法を指定できます。
• ジェネレータの組み合わせ: + や * 演算子を使って、複数のジェネレータを組み合わせることも可能です。例えば、定義域全体に一様に点をサンプリングしつつ、境界付近にはより多くの点を配置する、といった複雑なサンプリング戦略も実現できます。

さらに高度な機能（簡単な紹介）

neurodiffeq は、より進んだ応用や研究に対応するための機能も備えています。

• Solution Bundle (解の束) と Reverse Problems (逆問題):
  • Solution Bundle: 微分方程式に含まれるパラメータ（例えば、Lotka-Volterra方程式の \(\alpha\),\(\beta\),\(\gamma\),\(\delta\)）が異なる複数のケースについて、それぞれの解を個別に計算するのは非効率的な場合があります。BundleSolver を使うと、これらのパラメータをネットワークの追加の入力として扱い、パラメータが変化したときに解がどのように変わるかという「解の束（Solution Bundle）」を一度の学習で効率的に得ることができます。条件もパラメータ化された BundleIVP などを使用します。
  • Reverse Problems (Inverse Problems): 得られた Solution Bundle を利用して、観測データに最もよく適合する微分方程式のパラメータを推定する「逆問題」に取り組むことも可能です。これは、実験データから物理モデルの未知の定数を決定する際などに非常に有用です。
• Transfer Learning (転移学習): ある問題で学習済みのニューラルネットワークの重みを、別の関連する問題の初期値として利用する転移学習が可能です。これにより、新しい問題の学習時間を短縮したり、より良い解に到達したりすることが期待できます。neurodiffeq では、PyTorch の標準的なネットワーク保存・読み込み機能 (torch.save, torch.load) を利用して、学習済みネットワークを新しいソルバーに渡すことで転移学習を実現できます。

ハイパーパラメータ調整のヒント

ニューラルネットワークを用いた解法は強力ですが、その性能はハイパーパラメータの選択に大きく依存します。以下は、neurodiffeq を使う上でのいくつかのヒントです。

• ネットワーク構造: 隠れ層の数や各層のユニット数は、問題の複雑さに応じて調整が必要です。単純な問題であれば小さなネットワークで十分ですが、複雑な解を表現するにはより大きなネットワークが必要になることがあります。
• 活性化関数: neurodiffeq のドキュメントでは、ReLU 活性化関数は2階微分が0になってしまうため、高階微分方程式を解く際には避けるべきであると推奨されています。代わりに Tanh, SinActv, Swish (SiLU) などがよく用いられます。
• オプティマイザと学習率: PyTorch が提供する様々なオプティマイザ (Adam, SGDなど) を利用できます。学習率は学習の進捗に大きく影響するため、学習率スケジューラと合わせて調整することが一般的です。学習率の変更は、標準的なPyTorchの方法（オプティマイザの param_groups を操作する）で行います。
• サンプリング戦略: 前述のジェネレータを使い、損失を計算する点の分布や数を調整します。解が急激に変化する領域や、境界条件が重要な領域には、より多くの点をサンプリングすると良い結果が得られることがあります。
• 正規化・無次元化: 微分方程式や定義域を無次元化し、変数のスケールを 0 や 1 の付近に収めることは、学習の安定性や効率を向上させる上で非常に有効なテクニックです。
• カリキュラム学習: 非常に複雑な問題や広い定義域で解を求める場合、最初から全体を学習させようとするとうまくいかないことがあります。そのような場合は、まず狭い範囲や単純な条件で学習させ（簡単な問題）、徐々に範囲を広げたり条件を複雑にしたりする「カリキュラム学習」的なアプローチが有効な場合があります。

おわりに

本記事では、ニューラルネットワークを用いて微分方程式を解くためのPyTorchベースのライブラリ neurodiffeq について、その基本的な考え方から具体的な使用例、さらには高度な機能に至るまでを紹介しました。

neurodiffeq は、微分方程式という古典的かつ重要な問題に対し、深層学習という新しいアプローチを手軽に試すことができる強力なツールです。基本的な使い方をマスターすれば、ODEやPDEの数値解を連続関数として得られるだけでなく、Solution Bundleや逆問題といった応用的な課題にも挑戦できます。

もちろん、ニューラルネットワークに基づく解法は万能ではなく、ハイパーパラメータの調整や問題設定の工夫が求められる場面も多々あります。しかし、その柔軟性と表現力の高さは、従来の数値解法では困難だった問題への新たな道を開く可能性を秘めています。

今回の記事が、neurodiffeq を使った微分方程式への挑戦を始めるきっかけとなれば幸いです。ぜひ、ご自身の興味のある問題で、この新しい解法の可能性を探求してみてください。

More Information

• arXiv:2009.09943, Brandon Paulsen et al., 「NeuroDiff: Scalable Differential Verification of Neural Networks using Fine-Grained Approximation」, https://arxiv.org/abs/2009.09943
• arXiv:2201.03269, Andrea Sacchetti et al., 「Neural Networks to solve Partial Differential Equations: a Comparison with Finite Elements」, https://arxiv.org/abs/2201.03269
• arXiv:2405.20836, Chinmay Datar et al., 「Solving partial differential equations with sampled neural networks」, https://arxiv.org/abs/2405.20836
• arXiv:2502.12177, Shuheng Liu et al., 「Recent Advances of NeuroDiffEq — An Open-Source Library for Physics-Informed Neural Networks」, https://www.arxiv.org/abs/2502.12177