QUICK-PDE: Qiskit関数（ ColibriTD 著）

APIリファレンスを参照してください

概要

ここで紹介する偏微分方程式（PDE）ソルバーは、当社の量子革新的コンピューティングキット（QUICK）プラットフォーム（QUICK-PDE）の一部であり、Qiskit関数としてパッケージ化されています。 QUICK-PDE関数を使えば、IBM Quantum QPU上でドメイン固有の偏微分方程式を解くことができます。 この機能は、 ColibriTD's に記載されているH-DESの仕様書に記載されたアルゴリズムに基づいています。 このアルゴリズムは、計算流体力学（CFD）と材料変形（MD）を起点として、複雑なマルチフィジックス問題を解決でき、その他のユースケースも近日提供予定です。

微分方程式に取り組むために、試行解は直交関数（典型的にはチェビシェフ多項式、より具体的には $2^n$、 $n$ は関数をエンコードする量子ビットの数）の線形結合としてエンコードされ、可変量子回路（VQC）の角度によってパラメータ化される。 アナザッツは関数を符号化した状態を生成し、その関数はすべての点で関数を評価できるような組み合わせの観測量によって評価される。 次に、微分方程式がエンコードされた損失関数を評価し、次のようにハイブリッドループで角度を微調整することができる。 試行解は徐々に実際の解に近づいていき、満足のいく結果に到達する。

このハイブリッド・ループに加えて、異なるオプティマイザーを連鎖させることもできる。 これは、大域的なオプティマイザーで最適な角度のセットを見つけ、さらに微調整されたオプティマイザーで、隣接する角度の最適なセットへの勾配をたどりたい場合に便利です。 数値流体力学（CFD）の場合、デフォルトの最適化シーケンスで最良の結果が得られますが、材料変形（MD）の場合は、デフォルトで良好な結果が得られる一方で、問題固有の利点を得るためにさらに設定することができます。

この関数の各変数には、量子ビットの数を指定する（これは自由に弄ることができる）。 10個の同じ回路を積み重ね、1つの大きな回路を通して異なる量子ビットで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個のサンプルポイントがあります。

材料の変形

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

$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です。

ベンチマーク

以下の表は、我々の機能の様々な実行に関する統計である。

使用を開始する

QUICK-PDE 関数の利用申請を行うには、フォームに必要事項をご記入ください。 次に、 アカウントがすでにローカル環境に保存されていることを前提として、次のように関数を選択してください：

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)

サポートの利用

サポートについては、までご連絡ください [email protected]。

次のステップ