QUICK-PDE: Eine Qiskit Function von ColibriTD

Siehe die API-Referenz

hinweis

Qiskit Functions sind eine experimentelle Funktion, die Benutzern des IBM Quantum® Premium Plans, Flex Plans und On-Prem (über IBM Quantum Platform API) Plans zur Verfügung steht. Sie befinden sich im Preview-Release-Status und können sich ändern.

Überblick

Der hier vorgestellte Solver für partielle Differentialgleichungen (PDE) ist Teil unserer Quantum Innovative Computing Kit (QUICK)-Plattform (QUICK-PDE) und wird als Qiskit Function bereitgestellt. Mit der QUICK-PDE-Funktion kannst du domänenspezifische partielle Differentialgleichungen auf IBM Quantum QPUs lösen. Diese Funktion basiert auf dem in dem H-DES-Paper von ColibriTD beschriebenen Algorithmus. Dieser Algorithmus kann komplexe Multiphysik-Probleme lösen, beginnend mit Computational Fluid Dynamics (CFD) und Materialverformung (MD), wobei weitere Anwendungsfälle in Kürze folgen.

Um die Differentialgleichungen anzugehen, werden die Versuchslösungen als Linearkombinationen orthogonaler Funktionen (typischerweise Chebyshev-Polynome, und genauer $2^n$ davon, wobei $n$ die Anzahl der Qubits ist, die deine Funktion kodieren) kodiert, parametrisiert durch die Winkel eines variablen Quantum Circuits (Variable Quantum Circuit, VQC). Der Ansatz erzeugt einen Zustand, der die Funktion kodiert; dieser Zustand wird durch Observablen ausgewertet, deren Kombinationen die Auswertung der Funktion an allen Punkten ermöglichen. Du kannst dann die Verlustfunktion auswerten, in der die Differentialgleichungen kodiert sind, und die Winkel in einer hybriden Schleife feinabstimmen, wie im Folgenden gezeigt. Die Versuchslösungen kommen den tatsächlichen Lösungen schrittweise näher, bis du ein zufriedenstellendes Ergebnis erreichst.

Zusätzlich zu dieser hybriden Schleife kannst du auch verschiedene Optimierer verketten. Dies ist nützlich, wenn du einen globalen Optimierer verwenden möchtest, um einen guten Satz von Winkeln zu finden, und dann einen feineren Optimierer, um einem Gradienten zum besten Satz benachbarter Winkel zu folgen. Im Fall der numerischen Strömungsmechanik (CFD) liefert die Standard-Optimierungssequenz die besten Ergebnisse - aber bei der Materialverformung (MD) liefert die Standardsequenz zwar gute Ergebnisse, du kannst sie jedoch für problemspezifische Vorteile weiter anpassen.

Beachte, dass wir für jede Variable der Funktion die Anzahl der Qubits angeben (mit der du experimentieren kannst). Durch das Stapeln von 10 identischen Circuits und die Auswertung der 10 identischen Observablen auf verschiedenen Qubits über einen einzigen großen Circuit hinweg kannst du innerhalb des CMA-Optimierungsprozesses Rauschen mindern, indem du dich auf die Noise-Learner-Methode stützt, und die Anzahl der benötigten Shots erheblich reduzieren.

Numerische Strömungsmechanik

Die inviszide Burgers-Gleichung modelliert strömende nicht-viskose Flüssigkeiten wie folgt:

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

$u$ stellt das Fluidgeschwindigkeitsfeld dar. Dieser Anwendungsfall hat eine zeitliche Randbedingung: Du kannst die Anfangsbedingung auswählen und dann das System relaxieren lassen. Derzeit sind die einzigen akzeptierten Anfangsbedingungen lineare Funktionen: $ax + b$. Die analytische Lösung lautet:

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

Die druckfreien Euler-Gleichungen modellieren kompressible, nicht-viskose Fluidströmung mit Dämpfung wie folgt:

$\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$ stellt das Dichtefeld dar, $u$ das Geschwindigkeitsfeld und $\mu$ einen Dämpfungskoeffizienten. In unserer Formulierung setzen wir $\lambda = 1$, sodass es im Folgenden nicht als Parameter verwendet wird. Dieser Anwendungsfall hat zeitliche Randbedingungen: $g(0, x) = e^{-x}$ und $u(0, x) = x$. Die analytische Lösung lautet:

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

Die Argumente für CFDs Differentialgleichungen befinden sich auf einem festen Gitter wie folgt:

• $t$ liegt zwischen 0 und 0,95 mit 41 Abtastpunkten. $x$ liegt zwischen 0 und 0,95 mit 41 Abtastpunkten.

Materialverformung

Dieser Anwendungsfall konzentriert sich auf hypoelastische Verformung mit dem eindimensionalen Zugversuch, bei dem ein im Raum fixierter Stab an seinem anderen Ende gezogen wird. Wir beschreiben das Problem wie folgt:

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

$\sigma' - b = 0,$

$K$ stellt den Kompressionsmodul des gedehnten Materials dar, $n$ den Exponenten eines Potenzgesetzes, $b$ die Kraft pro Masseneinheit, $\epsilon_0$ die Proportionalitätsspannungsgrenze, $\sigma_0$ die Proportionalitätsdehnungsgrenze, $u$ die Spannungsfunktion und $\sigma$ die Dehnungsfunktion. Die analytische Lösung lautet:

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

wobei $\sigma_0 = g(0)$ die Randbedingung für die Dehnung bei $x=0$ ist.

Der betrachtete Stab hat eine Einheitslänge. Dieser Anwendungsfall hat eine Randbedingung für die Oberflächenlast (Randspannung) $t$, oder die Menge an Arbeit, die benötigt wird, um den Stab zu dehnen.

Die Argumente für MDs Differentialgleichungen befinden sich auf einem festen Gitter wie folgt:

• $x$ liegt zwischen 0 und 1 mit 30 Abtastpunkten.

Benchmarks

Die folgende Tabelle zeigt Statistiken zu verschiedenen Ausführungen unserer Funktion.

Beispiel	Anzahl der Qubits	Initialisierung	Fehler	Gesamtzeit (min)	Runtime-Nutzung (min)
Inviszide Burgers-Gleichung	50	PHYSICALLY_INFORMED	$10^{-2}$	66	25
Druckfreie Euler-Gleichungen	70	PHYSICALLY_INFORMED	$10^{-2}$	48	34
Hypoelastischer 1D-Zugversuch	18	RANDOM	$10^{-2}$	123	100

Erste Schritte

Fülle das Formular aus, um Zugang zur QUICK-PDE-Funktion anzufordern. Anschließend, vorausgesetzt du hast dein Konto bereits in deiner lokalen Umgebung gespeichert, wähle die Funktion wie folgt aus:

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

Beispiele

Um loszulegen, probiere eines der folgenden Beispiele:

Inviszide Burgers-Gleichung (CFD)

Mit der Anfangsbedingung $u(0,x) = x$ liefert die Burgers-Gleichung folgende Ergebnisse:

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

Prüfe den Status deiner Qiskit Function-Workload oder rufe die Ergebnisse wie folgt ab:

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

Druckfreie Euler-Gleichung (CFD)

Für die Euler-Gleichung, wenn Anfangsbedingungen auf $g(0, x) = e^{-x}$ und $u(0, x) = x$ gesetzt sind, für ein gegebenes $\mu$ (hier $\mu = 0.1$) und $\lambda = 1$, sind die Ergebnisse wie folgt:

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

Materialverformung

Der Anwendungsfall Materialverformung erfordert die physikalischen Parameter deines Materials und die angewandte Kraft wie folgt:

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

Das Folgende ist ein Beispiel dafür, wie du den Wert der Funktion für einen bestimmten Satz von Koordinaten erhältst:

# 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

Fehlermeldungen abrufen

Wenn der Status deiner Workload ERROR ist, verwende job.error_message(), um die Fehlermeldung zum Debuggen abzurufen, wie folgt:

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)

Support erhalten

Für Support kontaktiere qiskit-function-support@colibritd.com.

Nächste Schritte

Empfehlungen

• Fülle das Formular aus, um Zugang zur QUICK-PDE-Funktion anzufordern.
• Besuche die API-Referenz für diese Qiskit Function.
• Versuche, eine strömende nicht-viskose Flüssigkeit mit QUICK-PDE im Tutorial zu modellieren.
• Lies Jaffali, H., et al. (2025). H-DES: a Quantum-Classical Hybrid Differential Equation Solver. arXiv preprint arXiv:2410.01130.