QUICK-PDE: ColibriTDによるQiskit Function

備考

Qiskit Functionsは、IBM Quantum® Premium Plan、Flex Plan、およびOn-Prem（IBM Quantum Platform API経由）Planユーザーが利用できる実験的な機能です。プレビューリリースの状態であり、変更される可能性があります。

概要

ここで紹介する偏微分方程式（PDE）ソルバーは、Quantum Innovative Computing Kit（QUICK）プラットフォーム（QUICK-PDE）の一部であり、Qiskit Functionとしてパッケージ化されています。QUICK-PDE functionを使用すると、IBM Quantum QPU上でドメイン固有の偏微分方程式を解くことができます。このfunctionは、ColibriTDのH-DES説明論文に記載されているアルゴリズムに基づいています。このアルゴリズムは、計算流体力学（CFD）や材料変形（MD）を皮切りに、複雑なマルチフィジックス問題を解くことができ、他のユースケースも近日公開予定です。

微分方程式を解くために、試行解は直交関数の線形結合（通常はチェビシェフ多項式、より具体的には関数をエンコードするQubit数を $n$ とした場合の $2^n$ 個）としてエンコードされ、可変量子Circuit（VQC）の角度によってパラメータ化されます。ansatzは関数をエンコードする状態を生成し、その組み合わせによってすべての点で関数を評価できるオブザーバブルによって評価されます。そのあと、微分方程式がエンコードされた損失関数を評価し、以下に示すようなハイブリッドループで角度を微調整します。試行解は徐々に実際の解に近づき、満足のいく結果が得られるまで繰り返されます。

QUICK-PDE functionのワークフロー

このハイブリッドループに加えて、異なるオプティマイザーを連鎖させることもできます。これは、グローバルオプティマイザーで良い角度のセットを見つけ、その後、より細かく調整されたオプティマイザーで勾配に沿って最良の隣接角度のセットを探索したい場合に便利です。計算流体力学（CFD）の場合、デフォルトの最適化シーケンスが最良の結果をもたらします。一方、材料変形（MD）の場合、デフォルトでも良い結果が得られますが、問題に応じてさらに設定を調整し、性能を高めることもできます。

関数の各変数に対してQubit数を指定します（変更可能です）。10個の同一Circuitを重ねて、1つの大きなCircuit全体の異なるQubitで10個の同一オブザーバブルを評価することで、CMA最適化プロセス内でノイズ軽減を行い、ノイズ学習法を活用して必要なショット数を大幅に削減できます。

計算流体力学

非粘性バーガーズ方程式は、非粘性流体の流れを次のようにモデル化します：

$\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} = 0,$

$u$ は流体速度場を表します。このユースケースには時間的境界条件があります：初期条件を選択し、システムを緩和させることができます。現在、受け入れられる初期条件は線形関数のみです： $ax + b$ 。解析解は次のとおりです：

$u(t, x) = \frac{ax + b}{at + 1}.$

圧力なしオイラー方程式は、以下のように減衰を伴う圧縮性非粘性流体の流れをモデル化します。

$\frac{\partial g}{\partial t} + u\frac{\partial g}{\partial x} + g\frac{\partial u}{\partial x} = 0,$

$u\frac{\partial g}{\partial t} + g\frac{\partial u}{\partial t} + u^2\frac{\partial g}{\partial x} + 2gu\frac{\partial u}{\partial x} + \frac{\mu}{(1+t)^{\lambda}} g u = 0,$

$g$ は密度場、 $u$ は速度場、 $\mu$ は減衰係数を表します。この定式化では $\lambda = 1$ としているため、以降ではパラメータとして使用されません。この使用例には時間的境界条件があります: $g(0, x) = e^{-x}$ および $u(0, x) = x$ 。解析解は次のとおりです。

$g(t, x) = \frac{1-\mu}{(1+t)^{1-\mu} - \mu} \exp\!\left(\frac{(\mu-1)\, x}{(1+t)^{1-\mu} - \mu}\right),$

$u(t, x) = \frac{(1-\mu)\, x}{\left((1+t)^{1-\mu} - \mu\right)(1+t)^{\mu}}.$

CFDの微分方程式の引数は、次のように固定グリッド上にあります：

$t$ は0から0.95の間で41サンプル点。 $x$ は0から0.95の間で41サンプル点。

材料変形

このユースケースは、空間に固定された棒材がもう一方の端で引っ張られる一次元引張試験における亜弾性（hypoelastic）変形に焦点を当てています。問題は次のように記述されます：

$u' - \frac{\sigma}{3K} - \frac{2}{\sqrt{3}}\epsilon_0\left(\frac{\sigma'}{\sigma_0\sqrt{3}}\right)^n = 0,$

$\sigma' - b = 0,$

$K$ は引き伸ばされる材料の体積弾性率、 $n$ はべき乗則の指数、 $b$ は単位質量あたりの力、 $\epsilon_0$ は比例応力限界、 $\sigma_0$ は比例ひずみ限界、 $u$ は応力関数、 $\sigma$ はひずみ関数を表します。解析解は次のとおりです：

$\sigma(x) = \sigma_0 - bx,$

$u(x) = -\frac{3^{-(3+n)/2}}{2bK(1+n)\,\sigma_0^{n}}\Biggl[3^{(1+n)/2}b^2\sigma_0^n(1+n)x^2 - 2\cdot 3^{(1+n)/2}b\sigma_0^n(1+n)\sigma_0 x - 12\epsilon_0 K\sigma_0^{1+n}$ $- 12b\epsilon_0 K\sigma_0^n x\left(\frac{-bx+\sigma_0}{\sigma_0}\right)^n + 12\epsilon_0 K\sigma_0^{n+1}\left(\frac{-bx+\sigma_0}{\sigma_0}\right)^n - 12\epsilon_0 K \sigma_0^{1+n}\Biggr],$

ここで、 $\sigma_0 = g(0)$ は $x=0$ におけるひずみの境界条件です。

対象となる棒材は単位長さです。このユースケースには、表面応力 $t$ （棒材を引き伸ばすのに必要な仕事量）の境界条件があります。

MDの微分方程式の引数は、次のように固定グリッド上にあります：

$x$ は0から1の間で30サンプル点。

ベンチマーク

次の表は、このfunctionの様々な実行に関する統計を示しています。

例	Qubit数	初期化	エラー	合計時間（分）	ランタイム使用量（分）
非粘性バーガーズ方程式	50	`PHYSICALLY_INFORMED`	$10^{-2}$	66	25
圧力なしオイラー方程式	70	`PHYSICALLY_INFORMED`	$10^{-2}$	48	34
亜弾性一次元引張試験	18	`RANDOM`	$10^{-2}$	123	100

はじめに

QUICK-PDE functionへのアクセスリクエストフォームに記入してください。そのあと、すでにお使いのローカル環境にアカウントを保存済みであると仮定して、次のようにfunctionを選択します：

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"matplotlib": "matplotlib", "numpy": "numpy", "qiskit_ibm_catalog": "qiskit-ibm-catalog"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

from qiskit_ibm_catalog import QiskitFunctionsCatalog

catalog = QiskitFunctionsCatalog(
    channel="ibm_cloud / ibm_quantum_platform",
    instance="USER_CRN / HGP",
    token="USER_API_KEY / IQP_API_TOKEN",
)

catalog = QiskitFunctionsCatalog(channel="ibm_quantum_platform")

# Verify that you have access to the function
catalog.list()

quick = catalog.load("colibritd/quick-pde")

例

はじめに、次の例のいずれかを試してみてください：

非粘性バーガーズ方程式（CFD）

バーガーズ方程式について、初期条件を $u(0,x) = x$ に設定した場合、結果は次のようになります。

# launch the simulation with initial conditions u(0,x) = a*x + b
job = quick.run(
    use_case="CFD_BURGER", physical_parameters={"a": 1.0, "b": 0.0}
)

Qiskit Function ワークロードのステータスを確認するか、次のように結果を取得します：

# Print the ID so you can use it later, if necessary
print(job.job_id)
print(job.status())
solution = job.result()

import numpy as np
import matplotlib.pyplot as plt

def plot_result_3d(result):
    fig = plt.figure()
    ax = fig.add_subplot(projection="3d")

    t, x = np.meshgrid(result["samples"]["t"], result["samples"]["x"])

    ax.plot_surface(
        t,
        x,
        result["functions"]["u"],
        edgecolor="royalblue",
        lw=0.25,
        rstride=26,
        cstride=26,
        alpha=0.3,
    )
    ax.scatter(t, x, result["functions"]["u"], marker=".")
    ax.set(xlabel="t", ylabel="x", zlabel="u(t,x)")

    plt.show()

# Call
plot_result_3d(solution)

圧力なしオイラー方程式（CFD）

オイラー方程式について、初期条件を $g(0, x) = e^{-x}$ および $u(0, x) = x$ に設定し、与えられた $\mu$ （ここでは $\mu = 0.1$ ）と $\lambda = 1$ の場合、結果は次のようになります。

# Launches the solving for an arbitrary mu
job = quick.run(use_case="CFD_EULER", physical_parameters={"mu": 0.1})

solution = job.result()

# Colorplot function
def plot_result_2d(result):
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

    configs = {
        "g": {"cmap": "viridis", "title": "g(t, x)"},
        "u": {"cmap": "plasma", "title": "u(t, x)"},
    }

    t = result["samples"]["t"]
    x = result["samples"]["x"]

    for ax, (field, cfg) in zip(axes, configs.items()):
        v = result["functions"][field]

        im = ax.contourf(t, x, v, levels=50, cmap=cfg["cmap"])
        fig.colorbar(im, ax=ax, label=cfg["title"])

        ax.set_xlabel("t")
        ax.set_ylabel("x")
        ax.set_title(cfg["title"], fontsize=13, fontweight="bold")

    plt.tight_layout()
    plt.show()

plot_result_2d(solution)

材料変形

材料変形のユースケースでは、使用する材料の物理パラメータと印加する力が必要です。次のように指定します：

# Select the properties of your material
job = quick.run(
    use_case="MD",
    physical_parameters={
        "t": 12.0,
        "K": 100.0,
        "n": 4.0,
        "b": 10.0,
        "epsilon_0": 0.1,
        "sigma_0": 5.0,
    },
)

# Plot the result
solution = job.result()

_ = plt.figure()
stress_plot = plt.subplot(211)
plt.plot(solution["samples"]["x"], solution["functions"]["u"])
strain_plot = plt.subplot(212)
plt.plot(solution["samples"]["x"], solution["functions"]["sigma"])

plt.show()

特定の座標セットにおける関数の値を取得する方法の例を以下に示します：

# u(t=0.2, x=0.7) == 2
assert solution["samples"]["t"][1] == 0.2
assert solution["samples"]["x"][2] == 0.7
assert solution["functions"]["u"][1, 2] == 2

エラーメッセージの取得

ワークロードのステータスがERRORの場合、job.error_message()を使用してエラーメッセージを取得し、デバッグに役立ててください。次のように実行します：

job = quick.run(use_case="MD", physical_params={})

print(job.error_message())

# A wrapper can also be used for a more human readable version
def pprint_error(job):
    print("".join(eval(job.error_message())["error"]))

print("___")
pprint_error(job)

サポートを受ける

サポートについては、qiskit-function-support@colibritd.com にお問い合わせください。

次のステップ

推奨事項

QUICK-PDE functionへのアクセスリクエストフォームに記入してください。
このQiskit FunctionのAPIリファレンスを参照してください。
チュートリアルでQUICK-PDEを使用して非粘性流体の流れをモデル化してみてください。
Jaffali, H., et al. (2025). H-DES: a Quantum-Classical Hybrid Differential Equation Solver. arXiv preprint arXiv:2410.01130を参照してください。

概要​

計算流体力学​

材料変形​

ベンチマーク​

はじめに​

例​

非粘性バーガーズ方程式（CFD）​

圧力なしオイラー方程式（CFD）​

材料変形​

エラーメッセージの取得​

サポートを受ける​

次のステップ​

概要

計算流体力学

材料変形

ベンチマーク

はじめに

例

非粘性バーガーズ方程式（CFD）

圧力なしオイラー方程式（CFD）

材料変形

エラーメッセージの取得

サポートを受ける

次のステップ