QUICK-PDEを用いて流動する非粘性流体をモデル化する

注

Qiskit 関数は、 IBM Quantum® Premium Plan、Flex Plan、およびオンプレム ( IBM Quantum Platform API 経由) Plan ユーザーだけが利用できる実験的な機能です。これらはプレビューリリースであり、変更される可能性がある。

使用時間の目安：Heron r2 プロセッサーで50分。 (注：これはあくまでも目安です。ランタイムは異なるかもしれない)。

この関数の実行時間は通常20分以上かかることに注意してください、ですから、このチュートリアルを2つのセクションに分けた方がよいでしょう。このチュートリアルを2つのセクションに分けるとよいでしょう。 (ジョブが完了するのに十分な時間を与える)、ジョブの結果を使って作業する。

背景

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

マルチフィジックス問題（流体力学、熱拡散、材料変形など）は、偏流体力学、熱拡散、材料変形などのマルチフィジックス問題は、偏微分方程式(PDE)で記述することができる。微分方程式（PDE）で記述することができる。

このような問題は、さまざまな産業と密接に関連しており、応用数学の重要な一分野を構成している。応用数学の重要な一分野である。しかし多変量連立PDEを古典的なツールで解くことは、指数関数的に大きな資源を必要とするため、依然として困難である。指数関数的に大量のリソースを必要とするためである。

この関数は、複雑で変数が増える方程式に適しています。この関数は、複雑で変数が増える方程式に適しており、かつては難解とされていた最初のステップとなる。 PDEでモデル化された問題を完全に記述するには、初期条件と境界条件を知る必要がある。初期条件と境界条件を知る必要がある。これらはは、PDEの解とその解を求める経路を大きく変える。

このチュートリアルでは、その方法を紹介する：

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

要件

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

Qiskit SDK v2.0 またはそれ以降 (pip install qiskit)
Qiskit機能カタログ (pip install qiskit-ibm-catalog)
Matplotlib (pip install matplotlib)
QUICK-PDE 関数の呼び出し。アクセスをご希望の方は、フォームにご記入ください。

セットアップ

APIキーを使用して認証し、以下のように機能を選択します：

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

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

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

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

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

数値流体力学(CFD)の応用範囲は広い。そのため、基礎となるPDEを研究し解くことは重要である。重要なPDEsの一群ナビエ・ストークス方程式は、流体の運動を記述する非線形偏微分方程式系である。流体の運動を記述する非線形偏微分方程式系である。科学的問題や工学的応用に科学的な問題や工学的な応用に大きく関係している。

ある条件下では、ナビエ・ストークス方程式はブルガーの方程式に還元される。方程式に還元される。これは、流体力学、気体力学、非線形音響学などで起こる現象を記述する対流拡散方程式である。散逸系をモデル化する。

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

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

ここで、 $u(x,t)$ は、所定の位置 $x$ および時刻 $t$ における流体の流速場であり、 $\nu$ は流体の粘度である。粘度は流体の重要な性質の一つであり、その運動や変形に対する速度依存性の抵抗を表すものであり、したがって、流体の力学を決定する上で極めて重要な役割を果たしている。流体の粘性がゼロ（ $\nu = 0$ ）の場合、この方程式は保存方程式となり、その内部抵抗がないために不連続性（衝撃波）が生じることがある。この場合、この方程式は「非粘性バーガーズ方程式」と呼ばれ、非線形波動方程式の特例である。

厳密には、粘性流は自然界では発生しないが、空気力学的な流れをモデリングする場合空気力学的流れをモデル化する場合、輸送の影響が限りなく小さいため、問題の問題の非粘性記述は有用である。驚くことに、空力理論の70％以上は 70％以上の空力理論が非粘性流を扱っている。

このチュートリアルでは、CFDの例として非粘性バーガーズ方程式を用い、 QUICK-PDEを使用して IBM® のQPU上で以下の式を解きます：

$\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$ を調整して、それらが解法プロセスや解にどのような影響を与えるかを確認することができます。この初期条件を選択した場合、この方程式には以下の解析解が存在する：

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

これは、衝撃波が形成されていない限り成り立ちます。例えば、 $t < -1/a$ において、 $a < 0$ の場合などが挙げられます。

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

print(job.result())

Output:

{'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$ によって、変数 $t$ と $x$ ごとに必要な量子ビットが増えたり減ったりします。 $b$ .以下は、変数あたりの関数あたりの量子ビット数、変数あたりの関数あたりの深さ、変数あたりの関数あたりの量子ビット数を調整したものである。変数あたりの関数あたりの量子ビット数、関数あたりの深さ、ショット数を調整します。

バックエンドと実行モードの指定方法もご覧ください。

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

job_2 = quick.run(
    use_case="CFD_BURGER",
    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())

Output:

{'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: アルゴリズムの性能を比較する

当社のソリューション（HDES）の収束プロセスを、物理情報型ニューラルネットワーク（PINN）アルゴリズムおよびソルバー job_2 の性能と比較できます（論文および GitHub 関連リポジトリを参照）。

' job_2 s出力（量子ベースのアプローチ）の例では、古典ソルバーで最適化されるパラメータはわずか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個のパラメータの最適化を伴う。その後、ステップ315で目標損失が達成される。損失： 0.0014988397。

さて、公平な比較を行うためには、両方のケースで同じオプティマイザを使用すべきです。 12層×20ニューロンで確認した最小反復回数は、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

同様に、あなたのデータから同じ操作を行い、PINN解との job_2 比較プロットを作成できます。

# 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()

Output:

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()

Output:

チュートリアル調査

このチュートリアルに関するフィードバックをお寄せください。あなたの洞察は、私たちのコンテンツの提供とユーザーエクスペリエンスを向上させるのに役立ちます：

アンケートへのリンク

このページは役に立ちましたか？

バグや誤字の報告、またはコンテンツの要求はGitHubで行ってください。