QUICK-PDE: una Qiskit Function di ColibriTD

Consulta il riferimento API

nota

Le Qiskit Functions sono una funzionalità sperimentale disponibile per gli utenti dei piani IBM Quantum® Premium, Flex e On-Prem (tramite IBM Quantum Platform API). Sono in stato di anteprima e soggette a modifiche.

Panoramica

Il risolutore di Equazioni Differenziali alle Derivate Parziali (PDE) presentato qui fa parte della nostra piattaforma Quantum Innovative Computing Kit (QUICK) (QUICK-PDE) ed è distribuito come Qiskit Function. Con la funzione QUICK-PDE puoi risolvere equazioni differenziali alle derivate parziali specifiche per dominio sui QPU IBM Quantum. Questa funzione si basa sull'algoritmo descritto nel paper descrittivo H-DES di ColibriTD. Questo algoritmo è in grado di risolvere problemi multi-fisici complessi, a partire dalla Fluidodinamica Computazionale (CFD) e dalla Deformazione dei Materiali (MD), con altri casi d'uso in arrivo.

Per affrontare le equazioni differenziali, le soluzioni di tentativo sono codificate come combinazioni lineari di funzioni ortogonali (tipicamente polinomi di Chebyshev, e più specificamente $2^n$ di essi, dove $n$ è il numero di qubit che codificano la funzione), parametrizzate dagli angoli di un Variable Quantum Circuit (VQC). L'ansatz genera uno stato che codifica la funzione, valutata tramite osservabili le cui combinazioni consentono di calcolare la funzione in tutti i punti. Puoi quindi valutare la funzione di perdita in cui sono codificate le equazioni differenziali e affinare gli angoli in un ciclo ibrido, come mostrato di seguito. Le soluzioni di tentativo si avvicinano progressivamente alle soluzioni reali finché non si raggiunge un risultato soddisfacente.

Oltre a questo ciclo ibrido, puoi anche concatenare diversi ottimizzatori. Questo è utile quando vuoi che un ottimizzatore globale trovi un buon insieme di angoli, e poi un ottimizzatore più fine segua un gradiente verso il miglior insieme di angoli vicini. Nel caso della fluidodinamica computazionale (CFD), la sequenza di ottimizzazione predefinita produce i migliori risultati; nel caso della deformazione dei materiali (MD), il valore predefinito fornisce buoni risultati, ma puoi configurarlo ulteriormente per ottenere vantaggi specifici per il problema.

Nota che per ogni variabile della funzione specifichiamo il numero di qubit (su cui puoi sperimentare). Impilando 10 circuiti identici e valutando i 10 osservabili identici su qubit diversi all'interno di un unico grande circuito, puoi mitigare il rumore durante il processo di ottimizzazione CMA, sfruttando il metodo noise learner, e ridurre significativamente il numero di shot necessari.

Fluidodinamica computazionale

L'equazione di Burgers non viscosa modella fluidi non viscosi in moto nel modo seguente:

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

$u$ rappresenta il campo di velocità del fluido. Questo caso d'uso ha una condizione al contorno temporale: puoi selezionare la condizione iniziale e quindi lasciare che il sistema si rilasci. Attualmente le uniche condizioni iniziali accettate sono funzioni lineari: $ax + b$. La soluzione analitica è:

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

Le equazioni di Eulero senza pressione modellano il flusso di fluido inviscido comprimibile con smorzamento come segue:

$\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$ rappresenta il campo di densità, $u$ il campo di velocità e $\mu$ un coefficiente di smorzamento. Nella nostra formulazione, poniamo $\lambda = 1$, quindi non verrà usato come parametro in seguito. Questo caso d'uso ha condizioni al contorno temporali: $g(0, x) = e^{-x}$ e $u(0, x) = x$. La soluzione analitica è:

$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}}.$

Gli argomenti delle equazioni differenziali per CFD si trovano su una griglia fissa, come segue:

• $t$ è compreso tra 0 e 0.95 con 41 punti campione. $x$ è compreso tra 0 e 0.95 con 41 punti campione.

Deformazione dei materiali

Questo caso d'uso si concentra sulla deformazione ipoelastica con il test di trazione monodimensionale, in cui una barra fissata nello spazio viene tirata all'estremità opposta. Descriviamo il problema come segue:

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

$\sigma' - b = 0,$

$K$ rappresenta il modulo di compressione del materiale in trazione, $n$ l'esponente di una legge di potenza, $b$ la forza per unità di massa, $\epsilon_0$ il limite di proporzionalità della tensione, $\sigma_0$ il limite di proporzionalità della deformazione, $u$ la funzione di tensione e $\sigma$ la funzione di deformazione. La soluzione analitica è:

$\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],$

dove $\sigma_0 = g(0)$ è la condizione al contorno sulla deformazione a $x=0$.

La barra considerata ha lunghezza unitaria. Questo caso d'uso ha una condizione al contorno per la tensione superficiale $t$, ovvero la quantità di lavoro necessaria per allungare la barra.

Gli argomenti delle equazioni differenziali per MD si trovano su una griglia fissa, come segue:

• $x$ è compreso tra 0 e 1 con 30 punti campione.

Benchmark

La tabella seguente presenta statistiche su varie esecuzioni della nostra funzione.

Esempio	Numero di qubit	Inizializzazione	Errore	Tempo totale (min)	Utilizzo runtime (min)
Equazione di Burgers non viscosa	50	PHYSICALLY_INFORMED	$10^{-2}$	66	25
Equazioni di Eulero senza pressione	70	PHYSICALLY_INFORMED	$10^{-2}$	48	34
Test di trazione 1D ipoelastico	18	RANDOM	$10^{-2}$	123	100

Per iniziare

Compila il modulo per richiedere l'accesso alla funzione QUICK-PDE. Quindi, supponendo che tu abbia già salvato il tuo account nel tuo ambiente locale, seleziona la funzione come segue:

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

Esempi

Per iniziare, prova uno dei seguenti esempi:

Equazione di Burgers non viscosa (CFD)

Per l'equazione di Burgers, quando le condizioni iniziali sono impostate su $u(0,x) = x$, i risultati sono i seguenti:

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

Controlla lo stato del workload della tua Qiskit Function o recupera i risultati come segue:

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

Equazione di Eulero senza pressione (CFD)

Per l'equazione di Eulero, quando le condizioni iniziali sono impostate su $g(0, x) = e^{-x}$ e $u(0, x) = x$, per un dato $\mu$ (qui, $\mu = 0.1$) e $\lambda = 1$, i risultati sono i seguenti:

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

Deformazione dei materiali

Il caso d'uso della deformazione dei materiali richiede i parametri fisici del materiale e la forza applicata, come segue:

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

Il seguente è un esempio di come ottenere il valore della funzione per un insieme specifico di coordinate:

# 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

Recupero dei messaggi di errore

Se lo stato del workload è ERROR, usa job.error_message() per recuperare il messaggio di errore e facilitare il debug, come segue:

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)

Supporto

Per assistenza, contatta qiskit-function-support@colibritd.com.

Passi successivi

Raccomandazioni

• Compila il modulo per richiedere l'accesso alla funzione QUICK-PDE.
• Consulta il riferimento API per questa Qiskit Function.
• Prova a modellare un fluido non viscoso in moto con QUICK-PDE nel tutorial.
• Leggi Jaffali, H., et al. (2025). H-DES: a Quantum-Classical Hybrid Differential Equation Solver. arXiv preprint arXiv:2410.01130.