Skip to main content
IBM Quantum Platform

QUICK-PDE: Eine Qiskit-Funktion von ColibriTD

Siehe API-Referenz

Note

Qiskit-Funktionen sind eine experimentelle Funktion, die den Nutzern des IBM Quantum® Premium Plans, des Flex Plans und des On-Prem-Plans (über IBM Quantum Platform API) zur Verfügung steht. Sie befinden sich im Status einer Vorabversion und können sich noch ändern.


Übersicht

Der hier vorgestellte Solver für partielle Differentialgleichungen (PDE) ist Teil unserer Quantum Innovative Computing Kit (QUICK)-Plattform (QUICK-PDE) und als Qiskit-Funktion verfügbar. Mit der QUICK-PDE-Funktion können Sie domänenspezifische partielle Differentialgleichungen auf IBM Quantum QPUs lösen. Diese Funktion basiert auf dem Algorithmus, der in der Beschreibung von H-DES unter ColibriTD's beschrieben ist. Dieser Algorithmus kann komplexe multiphysikalische Probleme lösen, angefangen bei der numerischen Strömungsmechanik (CFD) und Materialverformung (MD) bis hin zu weiteren Anwendungsfällen, die in Kürze folgen werden.

Um die Differentialgleichungen zu lösen, werden die Versuchslösungen als Linearkombinationen orthogonaler Funktionen kodiert (typischerweise Tschebyscheff-Polynome, genauer gesagt 2n2^n, wobei nn die Anzahl der Qubits ist, die Ihre Funktion kodieren), die durch die Winkel eines variablen Quantenkreises (VQC) parametrisiert sind. Der Ansatz erzeugt einen Zustand, der die Funktion kodiert, die durch Beobachtungswerte ausgewertet wird, deren Kombinationen es ermöglichen, die Funktion in allen Punkten auszuwerten. Sie können dann die Verlustfunktion, in der die Differentialgleichungen kodiert sind, auswerten und die Winkel in einer hybriden Schleife fein abstimmen, wie im Folgenden gezeigt. Die Probelösungen nähern sich schrittweise den tatsächlichen Lösungen an, bis Sie ein zufriedenstellendes Ergebnis erzielen.

Arbeitsablauf der QUICK-PDE-Funktion

Zusätzlich zu dieser hybriden Schleife können Sie auch verschiedene Optimierer miteinander verketten. Dies ist nützlich, wenn Sie eine globale Optimierung wünschen, um einen guten Satz von Winkeln zu finden, und dann eine feiner abgestimmte Optimierung, um einem Gradienten zum besten Satz von benachbarten Winkeln zu folgen. Im Falle der numerischen Strömungsmechanik (CFD) liefert die Standardoptimierungssequenz die besten Ergebnisse. Im Falle der Materialverformung (MD) liefert die Standardsequenz zwar gute Ergebnisse, aber Sie können sie weiter konfigurieren, um problemspezifische Vorteile zu erzielen.

Beachten Sie, dass wir für jede Variable der Funktion die Anzahl der Qubits angeben (mit denen Sie spielen können). Durch das Stapeln von 10 identischen Schaltkreisen und die Auswertung der 10 identischen Observablen auf verschiedenen Qubits in einem großen Schaltkreis kann man im Rahmen des CMA-Optimierungsprozesses mit Hilfe der Noise-Learner-Methode das Rauschen reduzieren und die Anzahl der benötigten Schüsse erheblich verringern.

Numerische Strömungsmechanik

Die nichtviskose Burgers-Gleichung modelliert die Strömung nichtviskoser Fluide wie folgt:

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

uu stellt das Strömungsgeschwindigkeitsfeld dar. Dieser Anwendungsfall weist eine zeitliche Randbedingung auf: Sie können die Anfangsbedingung festlegen und das System anschließend sich selbst überlassen. Derzeit sind als Anfangsbedingungen ausschließlich lineare Funktionen zulässig: ax+bax + b. Die analytische Lösung lautet:

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

Die drucklosen Euler-Gleichungen modellieren die Strömung eines kompressiblen, inviskosen Fluids mit Dämpfung wie folgt:

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

u∂g∂t+g∂u∂t+u2∂g∂x+2gu∂u∂x+μ(1+t)λgu=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,

gg steht für das Dichtefeld, uu für das Geschwindigkeitsfeld und μ\mu für einen Dämpfungskoeffizienten. In unserer Formulierung setzen wir λ=1\lambda = 1 an, sodass dieser im Folgenden nicht als Parameter verwendet wird. Dieser Anwendungsfall weist zeitliche Randbedingungen auf: g(0,x)=e−xg(0, x) = e^{-x} und u(0,x)=xu(0, x) = x. Die analytische Lösung lautet:

g(t,x)=1−μ(1+t)1−μ−μexp⁡ ⁣((μ−1) x(1+t)1−μ−μ),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)=(1−μ) x((1+t)1−μ−μ)(1+t)μ.u(t, x) = \frac{(1-\mu)\, x}{\left((1+t)^{1-\mu} - \mu\right)(1+t)^{\mu}}.

Die Argumente für die CFD-Differentialgleichungen sind auf einem festen Gitter wie folgt:

  • tt liegt zwischen 0 und 0.95 mit 41 Messpunkten. xx liegt zwischen 0 und 0.95 mit 41 Messpunkten.

Materialverformung

Dieser Anwendungsfall konzentriert sich auf die hypoelastische Verformung im Rahmen eines eindimensionalen Zugversuchs, bei dem ein im Raum festgelegter Stab an seinem anderen Ende gezogen wird. Wir beschreiben das Problem wie folgt:

u′−σ3K−23ϵ0(σ′σ03)n=0,u' - \frac{\sigma}{3K} - \frac{2}{\sqrt{3}}\epsilon_0\left(\frac{\sigma'}{\sigma_0\sqrt{3}}\right)^n = 0,

σ′−b=0,\sigma' - b = 0,

KK steht für den Volumenmodul des gedehnten Materials, nn für den Exponenten eines Potenzgesetzes, bb für die Kraft pro Masseneinheit, ϵ0\epsilon_0 für die Grenze der proportionalen Spannung, σ0\sigma_0 für die Grenze der proportionalen Dehnung, uu für die Spannungsfunktion und σ\sigma für die Dehnungsfunktion. Die analytische Lösung lautet:

σ(x)=σ0−bx,\sigma(x) = \sigma_0 - bx,

u(x)=−3−(3+n)/22bK(1+n) σ0n[3(1+n)/2b2σ0n(1+n)x2−2⋅3(1+n)/2bσ0n(1+n)σ0x−12ϵ0Kσ01+nu(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ϵ0Kσ0nx(−bx+σ0σ0)n+12ϵ0Kσ0n+1(−bx+σ0σ0)n−12ϵ0Kσ01+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 σ0=g(0)\sigma_0 = g(0) die Randbedingung für die Dehnung am Punkt x=0x=0 ist.

Der betrachtete Stab ist von einheitlicher Länge. Für diesen Anwendungsfall gibt es eine Randbedingung für die Oberflächenspannung tt, d. h. den Arbeitsaufwand, der erforderlich ist, um den Stab zu dehnen.

Die Argumente für die Differentialgleichungen von MD sind auf einem festen Gitter wie folgt:

  • xx liegt zwischen 0 und 1 und umfasst 30 Messpunkte.

Benchmarks

Die folgende Tabelle enthält Statistiken zu verschiedenen Durchläufen unserer Funktion.

Beispiel
Anzahl der Qubits
Initialisierung
Fehler
Gesamtzeit (Min.)
Laufzeitverbrauch (min)
Unviskose Burgers'sche Gleichung50PHYSICALLY_INFORMED10−210^{-2}6625 GB
Drucklose Euler-Gleichungen70PHYSICALLY_INFORMED10−210^{-2}4834
Hypoelastischer 1D Zugversuch18.RANDOM10−210^{-2}123100

Erste Schritte

Füllen Sie das Formular aus, um Zugriff auf die QUICK-PDE-Funktion zu beantragen. Angenommen, Sie haben Ihr Konto bereits in Ihrer lokalen Umgebung gespeichert, wählen Sie die Funktion wie folgt aus:

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

Probieren Sie für den Anfang eines der folgenden Beispiele aus:

Die Gleichung von Inviscid Burgers (CFD)

Für die Burgers-Gleichung ergeben sich bei folgenden Anfangsbedingungen: u(0,x)=xu(0,x) = x 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}
)

Überprüfen Sie den Status Ihrer Qiskit-Funktions-Workload oder rufen Sie 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)

Drucklose Euler-Gleichung (CFD)

Für die Eulergleichung ergeben sich bei den Anfangsbedingungen g(0,x)=e−xg(0, x) = e^{-x} und u(0,x)=xu(0, x) = x sowie für die gegebenen Werte μ\mu (hier: μ=0.1\mu = 0.1 ) und λ=1\lambda = 1 folgende Ergebnisse:

# 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

Für den Anwendungsfall "Materialverformung" werden die physikalischen Parameter Ihres Materials und die angewandte Kraft wie folgt benötigt:

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

Im Folgenden finden Sie ein Beispiel dafür, wie Sie den Funktionswert für einen bestimmten Koordinatensatz ermitteln können:

# 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 Ihr Workload-Status ERROR lautet, verwenden Sie job.error_message() , um die Fehlermeldung zur Fehlersuche 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)

Unterstützung anfordern

Für Support wenden Sie sich bitte an [email protected].


Nächste Schritte

Empfehlungen
War diese Seite hilfreich?
Melden Sie auf GitHub einen Fehler, einen Tippfehler oder fordern Sie Inhalte an.