QUICK-PDE : une fonction Qiskit par ColibriTD

Consulte la référence API

remarque

Les Qiskit Functions sont une fonctionnalité expérimentale disponible pour les utilisateurs des plans IBM Quantum® Premium, Flex et On-Prem (via l'API IBM Quantum Platform). Elles sont en version préliminaire et susceptibles de changer.

Vue d'ensemble

Le solveur d'équations aux dérivées partielles (EDP) présenté ici fait partie de notre plateforme Quantum Innovative Computing Kit (QUICK) (QUICK-PDE), et est proposé sous forme de fonction Qiskit. Avec la fonction QUICK-PDE, tu peux résoudre des équations aux dérivées partielles spécifiques à un domaine sur des QPU IBM Quantum. Cette fonction est basée sur l'algorithme décrit dans l'article de ColibriTD décrivant H-DES. Cet algorithme peut résoudre des problèmes multi-physiques complexes, en commençant par la mécanique des fluides numérique (CFD) et la déformation des matériaux (MD), et d'autres cas d'utilisation à venir.

Pour traiter les équations différentielles, les solutions d'essai sont encodées comme des combinaisons linéaires de fonctions orthogonales (généralement des polynômes de Chebyshev, et plus spécifiquement $2^n$ d'entre eux où $n$ est le nombre de qubits encodant ta fonction), paramétrées par les angles d'un circuit quantique variable (VQC). L'ansatz génère un état encodant la fonction, qui est évalué par des observables dont les combinaisons permettent d'évaluer la fonction en tous les points. Tu peux ensuite évaluer la fonction de perte dans laquelle les équations différentielles sont encodées, et affiner les angles dans une boucle hybride, comme indiqué ci-après. Les solutions d'essai se rapprochent progressivement des solutions réelles jusqu'à ce que tu obtiennes un résultat satisfaisant.

En plus de cette boucle hybride, tu peux également enchaîner différents optimiseurs. C'est utile lorsque tu veux qu'un optimiseur global trouve un bon ensemble d'angles, puis un optimiseur plus fin pour suivre un gradient vers le meilleur ensemble d'angles voisins. Dans le cas de la mécanique des fluides numérique (CFD), la séquence d'optimisation par défaut produit les meilleurs résultats - alors que pour la déformation des matériaux (MD), si la séquence par défaut donne de bons résultats, tu peux encore l'affiner pour des avantages spécifiques au problème.

Note que pour chaque variable de la fonction, nous spécifions un nombre de qubits (que tu peux ajuster). En empilant 10 circuits identiques et en évaluant les 10 observables identiques sur différents qubits tout au long d'un grand circuit, tu peux atténuer le bruit dans le processus d'optimisation CMA, en t'appuyant sur la méthode d'apprentissage du bruit, et réduire significativement le nombre de shots nécessaires.

Mécanique des fluides numérique

L'équation de Burgers non visqueuse modélise les fluides non visqueux en écoulement comme suit :

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

$u$ représente le champ de vitesse du fluide. Ce cas d'utilisation a une condition aux limites temporelle : tu peux sélectionner la condition initiale puis laisser le système se relaxer. Actuellement, les seules conditions initiales acceptées sont les fonctions linéaires : $ax + b$. La solution analytique est :

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

Les équations d'Euler sans pression modélisent l'écoulement d'un fluide compressible non visqueux avec amortissement comme suit :

$\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$ représente le champ de densité, $u$ le champ de vitesse et $\mu$ un coefficient d'amortissement. Dans notre formulation, nous posons $\lambda = 1$, il ne sera donc pas utilisé comme paramètre par la suite. Ce cas d'usage a des conditions aux limites temporelles : $g(0, x) = e^{-x}$ et $u(0, x) = x$. La solution analytique est :

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

Les arguments pour les équations différentielles de la CFD sont sur une grille fixe, comme suit :

• $t$ est compris entre 0 et 0,95 avec 41 points d'échantillonnage. $x$ est compris entre 0 et 0,95 avec 41 points d'échantillonnage.

Déformation des matériaux

Ce cas d'utilisation se concentre sur la déformation hypoélastique avec le test de traction unidimensionnel, dans lequel une barre fixée dans l'espace est tirée à son autre extrémité. Nous décrivons le problème comme suit :

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

$\sigma' - b = 0,$

$K$ représente le module de compressibilité du matériau étiré, $n$ l'exposant d'une loi de puissance, $b$ la force par unité de masse, $\epsilon_0$ la limite de contrainte proportionnelle, $\sigma_0$ la limite de déformation proportionnelle, $u$ la fonction de contrainte et $\sigma$ la fonction de déformation. La solution analytique est :

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

où $\sigma_0 = g(0)$ est la condition aux limites sur la déformation à $x=0$.

La barre considérée est de longueur unitaire. Ce cas d'utilisation a une condition aux limites pour la contrainte de surface $t$, ou la quantité de travail nécessaire pour étirer la barre.

Les arguments pour les équations différentielles de la MD sont sur une grille fixe, comme suit :

• $x$ est compris entre 0 et 1 avec 30 points d'échantillonnage.

Benchmarks

Le tableau suivant présente des statistiques sur différentes exécutions de notre fonction.

Exemple	Nombre de qubits	Initialisation	Erreur	Temps total (min)	Utilisation runtime (min)
Équation de Burgers non visqueuse	50	PHYSICALLY_INFORMED	$10^{-2}$	66	25
Équations d'Euler sans pression	70	PHYSICALLY_INFORMED	$10^{-2}$	48	34
Test de traction hypoélastique 1D	18	RANDOM	$10^{-2}$	123	100

Premiers pas

Remplis le formulaire pour demander l'accès à la fonction QUICK-PDE. Ensuite, en supposant que tu as déjà enregistré ton compte dans ton environnement local, sélectionne la fonction comme suit :

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

Exemples

Pour démarrer, essaie l'un des exemples suivants :

Équation de Burgers non visqueuse (CFD)

Pour l'équation de Burgers, lorsque les conditions initiales sont définies à $u(0,x) = x$, les résultats sont les suivants :

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

Vérifie le statut de ta charge de travail Qiskit Function ou récupère les résultats comme suit :

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

Équation d'Euler sans pression (CFD)

Pour l'équation d'Euler, lorsque les conditions initiales sont définies à $g(0, x) = e^{-x}$ et $u(0, x) = x$, pour un $\mu$ donné (ici, $\mu = 0.1$) et $\lambda = 1$, les résultats sont les suivants :

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

Déformation des matériaux

Le cas d'utilisation de déformation des matériaux nécessite les paramètres physiques de ton matériau et la force appliquée, comme suit :

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

Voici un exemple montrant comment obtenir la valeur de la fonction pour un jeu de coordonnées donné :

# 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

Récupérer les messages d'erreur

Si le statut de ta charge de travail est ERROR, utilise job.error_message() pour récupérer le message d'erreur et aider au débogage, comme suit :

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)

Obtenir de l'aide

Pour obtenir de l'aide, contacte qiskit-function-support@colibritd.com.

Étapes suivantes

Recommandations

• Remplis le formulaire pour demander l'accès à la fonction QUICK-PDE.
• Consulte la référence API pour cette fonction Qiskit.
• Essaie de modéliser un fluide non visqueux en écoulement avec QUICK-PDE dans le tutoriel.
• Consulte Jaffali, H., et al. (2025). H-DES: a Quantum-Classical Hybrid Differential Equation Solver. arXiv preprint arXiv:2410.01130.