QUICK-PDEを使用して非粘性流体の流れをモデル化する

注意

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

使用量の目安: Heron r2プロセッサーで50分。（注: これはあくまで目安です。実際の実行時間は異なる場合があります。）

この関数の実行時間は一般的に20分以上かかるため、このチュートリアルを2つのセクションに分けて進めることをお勧めします。最初のセクションでは、チュートリアルを読みながらジョブを送信し、2つ目のセクションでは数時間後（ジョブが完了するのに十分な時間を確保した後）にジョブの結果を確認します。

背景

このチュートリアルでは、ColibriTDのH-DES（ハイブリッド微分方程式ソルバー）を使用して、156Q Heron R2 QPU上で複雑なマルチフィジックス問題を解くためのQUICK-PDE関数の使い方を入門レベルで説明します。 基盤となるアルゴリズムはH-DESの論文に記載されています。 なお、このソルバーは非線形方程式も解くことができます。

流体力学、熱拡散、材料変形など、マルチフィジックス問題はあらゆる場面で偏微分方程式（PDE）によって記述されます。

このような問題は、さまざまな産業において非常に重要であり、応用数学の重要な分野を構成しています。しかし、非線形の多変数連立偏微分方程式を古典的なツールで解くことは、指数関数的に大きなリソースが必要となるため、依然として困難です。

この関数は、複雑さと変数が増大する方程式に適しており、かつてはコンピュータでは計算不可能と考えられていた可能性を解き放つ第一歩です。偏微分方程式でモデル化された問題を完全に記述するには、初期条件と境界条件を知る必要があります。これらの条件は、PDEの解および解を求めるまでの経路を大きく変える可能性があります。

このチュートリアルでは、以下の内容を学びます:

1. 初期条件関数のパラメータを定義する。
2. 量子ビット数（微分方程式の関数をエンコードするために使用）、回路の深さ、およびショット数を調整する。
3. QUICK-PDEを実行して基盤となる微分方程式を解く。

要件

このチュートリアルを始める前に、以下がインストールされていることを確認してください:

• Qiskit SDK v2.0以降（pip install qiskit）
• Qiskit Functions Catalog（pip install qiskit-ibm-catalog）
• Matplotlib（pip install matplotlib）
• QUICK-PDE関数へのアクセス。アクセスリクエストフォームに記入してください。

セットアップ

APIキーを使用して認証を行い、以下のように関数を選択します:

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit-ibm-catalog

import numpy as np
import matplotlib.pyplot as plt
from qiskit_ibm_catalog import QiskitFunctionsCatalog

catalog = QiskitFunctionsCatalog(
    channel="ibm_quantum_platform",
    instance="INSTANCE_CRN",
    token="YOUR_API_KEY",  # Use the 44-character API_KEY you created and saved from the IBM Quantum Platform Home dashboard
)

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

ステップ1: 解く問題のプロパティを設定する

このチュートリアルでは、初期条件によって決定される物理問題と、量子コンピュータ上で流体力学の例を解くためのアルゴリズム的側面という2つの観点からユーザー体験を扱います。

数値流体力学（CFD）は幅広い応用範囲を持つため、基盤となるPDEを研究し解くことが重要です。PDEの重要なファミリーの1つがナビエ・ストークス方程式であり、これは流体の運動を記述する非線形偏微分方程式の連立系です。科学的問題や工学的応用において非常に重要です。

特定の条件下では、ナビエ・ストークス方程式はバーガース方程式に帰着されます。これは、散逸系をモデル化することで流体力学、気体力学、非線形音響学などで生じる現象を記述する移流拡散方程式です。

この方程式の1次元バージョンは、2つの変数に依存します: $t \in \mathbb{R}_{\geq 0}$は時間次元をモデル化し、$x \in \mathbb{R}$は空間次元を表します。方程式の一般的な形式は粘性バーガース方程式と呼ばれ、次のように表されます:

$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = \nu \frac{\partial^2 u}{\partial^2 x},$

ここで、$u(x,t)$は与えられた位置$x$と時刻$t$における流体の速度場であり、$\nu$は流体の粘性です。粘性は流体の重要な性質であり、運動や変形に対する速度依存の抵抗を測定するものであり、したがって流体の挙動の決定に重要な役割を果たします。流体の粘性がゼロの場合（$\nu$ = 0）、方程式は保存則となり、内部抵抗がないために不連続面（衝撃波）が発生する可能性があります。この場合、方程式は非粘性バーガース方程式と呼ばれ、非線形波動方程式の特殊なケースです。

厳密に言えば、非粘性流れは自然界では発生しませんが、空気力学的な流れをモデル化する際、輸送の影響が無限小であることから、問題を非粘性として記述することが有用な場合があります。驚くべきことに、空気力学理論の70%以上が非粘性流れを扱っています。

このチュートリアルでは、QUICK-PDEを使用してIBM® QPU上で解くCFDの例として、以下の式で表される非粘性バーガース方程式を使用します:

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

この問題の初期条件は線形関数に設定されます: $u(t=0,x) = ax + b,\text{ with }a,b\in\mathbb{R}$ ここで、$a$と$b$は解の形状に影響を与える任意の定数です。$a$と$b$を調整して、それらが解法プロセスと解にどのように影響するかを確認できます。

job = quick.run(
    use_case="cfd",
    physical_parameters={"a": 1.0, "b": 1.0},
)

print(job.result())

{'functions': {'u': array([[1.        , 0.96112378, 0.9230742 , 0.88616096, 0.85058445,
        0.81644741, 0.78376878, 0.75249908, 0.72253689, 0.69374562,
        0.66597013, 0.63905258, 0.61284684, 0.58723093, 0.56211691,
        0.53745752, 0.51324915, 0.48953036, 0.46637547, 0.44388257,
        0.4221554 , 0.40127848, 0.38128488, 0.36211604, 0.34357308,
        0.32525895, 0.30651089, 0.28632252, 0.26325504, 0.23533692],
       [1.2375    , 1.19267729, 1.14850734, 1.10544526, 1.06382155,
        1.02385326, 0.98565757, 0.94926734, 0.91464784, 0.88171402,
        0.85034771, 0.82041411, 0.79177677, 0.76431068, 0.73791248,
        0.71250742, 0.68805224, 0.66453346, 0.64196021, 0.62035121,
        0.59971506, 0.5800232 , 0.56117499, 0.54295419, 0.52497612,
        0.50662498, 0.48698059, 0.4647339 , 0.43809065, 0.40466247],
       [1.475     , 1.4242308 , 1.37394048, 1.32472956, 1.27705866,
        1.23125911, 1.18754636, 1.1460356 , 1.10675879, 1.06968242,
        1.03472529, 1.00177563, 0.9707067 , 0.94139043, 0.91370806,
        0.88755732, 0.86285533, 0.83953655, 0.81754494, 0.79681986,
        0.77727473, 0.75876792, 0.74106511, 0.72379234, 0.70637915,
        0.687991  , 0.66745028, 0.64314527, 0.61292625, 0.57398802],
       [1.7125    , 1.65578431, 1.59937362, 1.54401386, 1.49029576,
        1.43866495, 1.38943515, 1.34280386, 1.29886974, 1.25765082,
        1.21910288, 1.18313715, 1.14963664, 1.11847019, 1.08950364,
        1.06260722, 1.03765842, 1.01453964, 0.99312968, 0.97328851,
        0.95483439, 0.93751264, 0.92095522, 0.90463049, 0.88778219,
        0.86935702, 0.84791997, 0.82155665, 0.78776186, 0.74331358],
       [1.95      , 1.88733782, 1.82480676, 1.76329816, 1.70353287,
        1.6460708 , 1.59132394, 1.53957212, 1.49098069, 1.44561922,
        1.40348046, 1.36449867, 1.32856657, 1.29554994, 1.26529921,
        1.23765712, 1.21246152, 1.18954273, 1.16871442, 1.14975716,
        1.13239406, 1.11625736, 1.10084533, 1.08546864, 1.06918523,
        1.05072304, 1.02838966, 0.99996803, 0.96259746, 0.91263913]])}, 'samples': {'t': array([0.        , 0.03275862, 0.06551724, 0.09827586, 0.13103448,
       0.1637931 , 0.19655172, 0.22931034, 0.26206897, 0.29482759,
       0.32758621, 0.36034483, 0.39310345, 0.42586207, 0.45862069,
       0.49137931, 0.52413793, 0.55689655, 0.58965517, 0.62241379,
       0.65517241, 0.68793103, 0.72068966, 0.75344828, 0.7862069 ,
       0.81896552, 0.85172414, 0.88448276, 0.91724138, 0.95      ]), 'x': array([0.    , 0.2375, 0.475 , 0.7125, 0.95  ])}}

ステップ2（必要に応じて）: 量子ハードウェア実行のために問題を最適化する

デフォルトでは、ソルバーは物理情報に基づくパラメータを使用します。これは、与えられた量子ビット数と深さに対する初期回路パラメータであり、ソルバーはここから最適化を開始します。

ショット数もデフォルト値を持つパラメータの一部であり、微調整が重要です。

解こうとしている設定によっては、満足のいく解を得るためにアルゴリズムのパラメータを調整する必要がある場合があります。例えば、$a$と$b$に応じて、変数$t$および$x$ごとに必要な量子ビット数が増減する可能性があります。以下では、関数ごとの変数あたりの量子ビット数、関数ごとの深さ、およびショット数を調整します。

バックエンドと実行モードの指定方法も確認できます。

また、物理情報に基づくパラメータが最適化プロセスを誤った方向に導く場合があります。その場合、initialization戦略を"RANDOM"に設定することで無効にできます。

job_2 = quick.run(
    use_case="cfd",
    physical_parameters={"a": 0.5, "b": 0.25},
    nb_qubits={"u": {"t": 2, "x": 1}},
    depth={"u": 3},
    shots=[500, 2500, 5000, 10000],
    initialization="RANDOM",
    backend="ibm_kingston",
    mode="session",
)

print(job_2.result())

{'functions': {'u': array([[0.25      , 0.24856543, 0.24687708, 0.2449444 , 0.24277686,
        0.24038389, 0.23777496, 0.23495952, 0.23194702, 0.22874691,
        0.22536866, 0.22182171, 0.21811551, 0.21425952, 0.2102632 ,
        0.20613599, 0.20188736, 0.19752675, 0.19306361, 0.18850741,
        0.18386759, 0.1791536 , 0.17437491, 0.16954096, 0.16466122,
        0.15974512, 0.15480213, 0.1498417 , 0.14487328, 0.13990632],
       [0.36875   , 0.36681313, 0.36457201, 0.36203594, 0.35921422,
        0.35611615, 0.35275103, 0.34912817, 0.34525687, 0.34114643,
        0.33680614, 0.33224532, 0.32747327, 0.32249928, 0.31733266,
        0.31198271, 0.30645873, 0.30077002, 0.29492589, 0.28893564,
        0.28280857, 0.27655397, 0.27018116, 0.26369944, 0.2571181 ,
        0.25044645, 0.24369378, 0.23686941, 0.22998264, 0.22304275],
       [0.4875    , 0.48506084, 0.48226695, 0.47912748, 0.47565158,
        0.47184841, 0.46772711, 0.46329683, 0.45856672, 0.45354594,
        0.44824363, 0.44266894, 0.43683103, 0.43073904, 0.42440212,
        0.41782942, 0.4110301 , 0.4040133 , 0.39678818, 0.38936388,
        0.38174955, 0.37395435, 0.36598742, 0.35785791, 0.34957498,
        0.34114777, 0.33258544, 0.32389713, 0.315092  , 0.30617919],
       [0.60625   , 0.60330854, 0.59996188, 0.59621902, 0.59208895,
        0.58758067, 0.58270318, 0.57746549, 0.57187658, 0.56594545,
        0.55968112, 0.55309256, 0.54618879, 0.53897879, 0.53147158,
        0.52367614, 0.51560147, 0.50725658, 0.49865046, 0.48979211,
        0.48069053, 0.47135472, 0.46179367, 0.45201638, 0.44203186,
        0.4318491 , 0.42147709, 0.41092485, 0.40020136, 0.38931562],
       [0.725     , 0.72155625, 0.71765682, 0.71331056, 0.70852631,
        0.70331293, 0.69767926, 0.69163414, 0.68518643, 0.67834497,
        0.6711186 , 0.66351618, 0.65554655, 0.64721855, 0.63854104,
        0.62952285, 0.62017284, 0.61049986, 0.60051274, 0.59022035,
        0.57963151, 0.56875509, 0.55759992, 0.54617486, 0.53448874,
        0.52255042, 0.51036875, 0.49795257, 0.48531072, 0.47245205]])}, 'samples': {'t': array([0.        , 0.03275862, 0.06551724, 0.09827586, 0.13103448,
       0.1637931 , 0.19655172, 0.22931034, 0.26206897, 0.29482759,
       0.32758621, 0.36034483, 0.39310345, 0.42586207, 0.45862069,
       0.49137931, 0.52413793, 0.55689655, 0.58965517, 0.62241379,
       0.65517241, 0.68793103, 0.72068966, 0.75344828, 0.7862069 ,
       0.81896552, 0.85172414, 0.88448276, 0.91724138, 0.95      ]), 'x': array([0.    , 0.2375, 0.475 , 0.7125, 0.95  ])}}

ステップ3: アルゴリズムの性能を比較する

job_2の解（HDES）の収束過程を、物理情報に基づくニューラルネットワーク（PINN）アルゴリズムおよびソルバーの性能と比較できます（論文および関連するGitHubリポジトリを参照してください）。

job_2の出力（量子ベースのアプローチ）の例では、古典的なソルバーで最適化されるパラメータはわずか13個（12個の回路パラメータと1個のスケーリングパラメータ）です。 収束過程は以下の通りです:

optimizers:
   CMA: {'ftarget': np.float64(0.1), 'verb_disp': 10, 'maxiter': 100}
   CMA: {'ftarget': np.float64(0.005), 'verb_disp': 10, 'maxiter': 20}
   CMA: {'ftarget': np.float64(0.0025), 'verb_disp': 10, 'maxiter': 30}
   CMA: {'ftarget': np.float64(0.0005), 'verb_disp': 10, 'maxiter': 10}

500 shots
================== CMA =================
option:  {'ftarget': np.float64(0.1), 'verb_disp': 10, 'maxiter': 100}
0/100, loss: 0.02456641

1000 shots
================== CMA =================
option:  {'ftarget': np.float64(0.005), 'verb_disp': 10, 'maxiter': 20}
0/20, loss: 0.03641833
1/20, loss: 0.02461719
2/20, loss: 0.0283689
3/20, loss: 0.009898383
4/20, loss: 0.04454522
5/20, loss: 0.007019971
6/20, loss: 0.00811147
7/20, loss: 0.01592619
8/20, loss: 0.00764708
9/20, loss: 0.01401516
10/20, loss: 0.01767467
11/20, loss: 0.01220387

5000 shots
================== CMA =================
option:  {'ftarget': np.float64(0.0025), 'verb_disp': 10, 'maxiter': 30}
0/30, loss: 0.01024792
1/30, loss: 0.004343748
2/30, loss: 0.01450951
3/30, loss: 0.008591284
4/30, loss: 0.00266414
5/30, loss: 0.007923613
6/30, loss: 0.02023853
7/30, loss: 0.01031438
8/30, loss: 0.009513116
9/30, loss: 0.008132266
10/30, loss: 0.005787766
11/30, loss: 0.00390582

10000 shots
================== CMA =================
option:  {'ftarget': np.float64(0.0005), 'verb_disp': 10, 'maxiter': 10}
0/10, loss: 0.002386168
1/10, loss: 0.004024823
2/10, loss: 0.001311999
3/10, loss: 0.003433991
4/10, loss: 0.002339664
5/10, loss: 0.002978438
6/10, loss: 0.005458391
7/10, loss: 0.002026701
8/10, loss: 0.00207467
9/10, loss: 0.001947627
final_loss: 0.00151994463476429

これは、28回の反復の後に0.0015未満の損失に到達でき、わずかな古典パラメータの最適化のみで実現できることを意味します。

次に、論文で提案されたデフォルト設定で勾配ベースのオプティマイザーを使用するPINNの解と同じ比較を行います。最適化対象の13個のパラメータを持つ回路に相当するのがニューラルネットワークであり、少なくとも20ニューロンの8層が必要で、3021個のパラメータの最適化が必要です。 その場合、ターゲット損失はStep 315で到達し、loss: 0.0014988397となります。

ここで、公正な比較を行うために、両方のケースで同じオプティマイザーを使用する必要があります。 20ニューロンの12層 = 4701パラメータで見つかった最小の反復回数は以下の通りです:

(10_w,20)-aCMA-ES (mu_w=5.9,w_1=27%) in dimension 4701 (seed=351961)
Iterat #Fevals   function value  axis ratio  sigma  min&max std  t[m:s]
    1     20 5.398521572351456e-02 1.0e+00 9.98e-03  1e-02  1e-02 0:02.3
    2     40 5.444650724530220e-02 1.0e+00 9.97e-03  1e-02  1e-02 0:05.1
    3     60 4.447407275438309e-02 1.0e+00 9.95e-03  1e-02  1e-02 0:08.2
    4     80 2.068969979882240e-02 1.0e+00 9.94e-03  1e-02  1e-02 0:11.7
    6    120 1.028892211616039e-02 1.0e+00 9.91e-03  1e-02  1e-02 0:20.1
    7    140 5.140972323715687e-03 1.0e+00 9.90e-03  1e-02  1e-02 0:25.4
    9    180 3.811701666563749e-03 1.0e+00 9.87e-03  1e-02  1e-02 0:37.4
   10    200 3.189878538250923e-03 1.0e+00 9.85e-03  1e-02  1e-02 0:44.2
   12    240 2.547040116041899e-03 1.0e+00 9.83e-03  1e-02  1e-02 0:59.7
   14    280 2.166548743844032e-03 1.0e+00 9.80e-03  1e-02  1e-02 1:18.0
   15    300 1.783065614290535e-03 1.0e+00 9.79e-03  1e-02  1e-02 1:28.4
   16    320 2.045844215899706e-03 1.0e+00 9.78e-03  1e-02  1e-02 1:39.8
Stopping early: loss 0.001405 <= target 0.0015
CMA-ES finished. Best loss: 0.001404788694344461

job_2のデータでも同様の処理を行い、PINNの解との比較をプロットできます。

# check the loss function and compare between the two approaches
print(job_2.logs())

ステップ4: 結果を利用する

得られた解を使って、さまざまな処理を行うことができます。以下では、結果をプロットする方法を示します。

solution = job.result()

# Plot the solution of the second simulation job_2
_ = plt.figure()
ax = plt.axes(projection="3d")

# plot the solution using the 3d plotting capabilities of pyplot
t, x = np.meshgrid(solution["samples"]["t"], solution["samples"]["x"])
ax.plot_surface(
    t,
    x,
    solution["functions"]["u"],
    edgecolor="royalblue",
    lw=0.25,
    rstride=26,
    cstride=26,
    alpha=0.3,
)
ax.scatter(t, x, solution, marker=".")
ax.set(xlabel="t", ylabel="x", zlabel="u(t,x)")

plt.show()

2回目の実行における初期条件の違いと、それが結果に与える影響に注目してください:

solution_2 = job_2.result()

# Plot the solution of the second simulation job_2
_ = plt.figure()
ax = plt.axes(projection="3d")

# plot the solution using the 3d plotting capabilities of pyplot
t, x = np.meshgrid(solution_2["samples"]["t"], solution_2["samples"]["x"])
ax.plot_surface(
    t,
    x,
    solution_2["functions"]["u"],
    edgecolor="royalblue",
    lw=0.25,
    rstride=26,
    cstride=26,
    alpha=0.3,
)
ax.scatter(t, x, solution_2, marker=".")
ax.set(xlabel="t", ylabel="x", zlabel="u(t,x)")

plt.show()

チュートリアルアンケート

このチュートリアルに関するフィードバックをお寄せください。皆様のご意見は、コンテンツとユーザー体験の改善に役立てさせていただきます:

アンケートへのリンク

Note: This survey is provided by IBM Quantum and relates to the original English content. To give feedback on doQumentation's website, translations, or code execution, please open a GitHub issue.