QUICK-PDE: функція Qiskit від ColibriTD

Дивись довідник API

примітка

Qiskit Functions — це експериментальна можливість, доступна користувачам IBM Quantum® Premium Plan, Flex Plan та On-Prem (через IBM Quantum Platform API) Plan. Вони перебувають у статусі попереднього релізу та можуть змінюватись.

Огляд

Розв'язувач диференціальних рівнянь у частинних похідних (PDE), представлений тут, є частиною нашої платформи Quantum Innovative Computing Kit (QUICK) (QUICK-PDE) і запакований як функція Qiskit. За допомогою функції QUICK-PDE ти можеш розв'язувати предметно-орієнтовані диференціальні рівняння в частинних похідних на QPU IBM Quantum. Ця функція базується на алгоритмі, описаному в статті ColibriTD про H-DES. Цей алгоритм може розв'язувати складні мультифізичні задачі, починаючи з обчислювальної гідродинаміки (CFD) та деформації матеріалів (MD), і незабаром з'являться нові варіанти використання.

Щоб впоратися з диференціальними рівняннями, пробні розв'язки кодуються як лінійні комбінації ортогональних функцій (зазвичай поліноми Чебишева, і конкретніше $2^n$ таких, де $n$ — кількість кубітів, що кодують твою функцію), параметризовані кутами варіаційної квантової схеми (VQC). Анзац генерує стан, що кодує функцію, яку обчислюють спостережувані величини, комбінації яких дозволяють обчислити функцію у всіх точках. Потім ти можеш обчислити функцію втрат, у якій закодовано диференціальні рівняння, та налаштувати кути в гібридному циклі, як показано нижче. Пробні розв'язки поступово наближаються до фактичних розв'язків, доки не буде досягнуто задовільного результату.

На додаток до цього гібридного циклу, ти також можеш поєднувати різні оптимізатори в ланцюжок. Це корисно, коли потрібно, щоб глобальний оптимізатор знайшов хороший набір кутів, а потім точніший оптимізатор пішов за градієнтом до найкращого набору сусідніх кутів. У випадку обчислювальної гідродинаміки (CFD) стандартна послідовність оптимізації дає найкращі результати — але у випадку деформації матеріалів (MD), хоча стандартні налаштування і дають хороші результати, ти можеш налаштувати її додатково для переваг, специфічних для конкретної задачі.

Зауваж, що для кожної змінної функції ми вказуємо кількість кубітів (з якою ти можеш експериментувати). Стекуючи 10 ідентичних схем та обчислюючи 10 ідентичних спостережуваних на різних кубітах в одній великій схемі, ти можеш знижувати шуми в процесі CMA-оптимізації, спираючись на метод навчання шумів, та значно зменшити потрібну кількість знімків (shots).

Обчислювальна гідродинаміка

Невʼязке рівняння Бюргерса моделює потік невʼязкої рідини таким чином:

$\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 точками вибірки.

Бенчмарки

У наступній таблиці представлено статистику різних запусків нашої функції.

Приклад	Кількість кубітів	Ініціалізація	Похибка	Загальний час (хв)	Використання Runtime (хв)
Inviscid Burgers' equation	50	PHYSICALLY_INFORMED	$10^{-2}$	66	25
Pressureless Euler's equations	70	PHYSICALLY_INFORMED	$10^{-2}$	48	34
Hypoelastic 1D tensile test	18	RANDOM	$10^{-2}$	123	100

Початок роботи

Заповни форму для запиту доступу до функції QUICK-PDE. Потім, якщо ти вже зберіг свій обліковий запис у локальному середовищі, вибери функцію таким чином:

# 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 таким чином:

# 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.
• Відвідай довідник API для цієї функції Qiskit.
• Спробуй змоделювати потік невʼязкої рідини за допомогою QUICK-PDE у підручнику.
• Ознайомся з Jaffali, H., та ін. (2025). H-DES: a Quantum-Classical Hybrid Differential Equation Solver. arXiv preprint arXiv:2410.01130.