QUICK-PDE : une fonction Qiskit par ColibriTD

remarque

Les fonctions Qiskit 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 sont susceptibles d'évoluer.

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 packagé 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 les QPU IBM Quantum. Cette fonction repose sur l'algorithme décrit dans l'article de présentation H-DES de ColibriTD. Cet algorithme peut résoudre des problèmes multi-physiques complexes, en commençant par la Dynamique des Fluides Numérique (CFD) et la Déformation des Matériaux (MD), avec d'autres cas d'usage à 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 précisément $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 illustré ci-dessous. Les solutions d'essai se rapprochent progressivement des solutions réelles jusqu'à obtenir 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 qu'un optimiseur plus précis suive un gradient vers le meilleur ensemble d'angles voisins. Dans le cas de la dynamique des fluides numérique (CFD), la séquence d'optimisation par défaut produit les meilleurs résultats — mais dans le cas de la déformation des matériaux (MD), bien que les valeurs par défaut donnent de bons résultats, tu peux la configurer davantage pour des bénéfices propres au problème.

À noter que pour chaque variable de la fonction, on spécifie un nombre de qubits (que tu peux ajuster). En empilant 10 circuits identiques et en évaluant les 10 observables identiques sur des qubits différents dans un seul grand circuit, tu peux atténuer le bruit dans le cadre du processus d'optimisation CMA, en s'appuyant sur la méthode de l'apprentissage du bruit, et réduire significativement le nombre de shots nécessaires.

Entrée/sortie

Dynamique des Fluides Numérique

L'équation de Burgers non visqueuse modélise l'écoulement de fluides non visqueux de la façon suivante :

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

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

Les arguments des équations différentielles de CFD sont définis sur une grille fixe, comme suit :

• $t$ est compris entre 0 et 0,95 avec 30 points d'échantillonnage. $x$ est compris entre 0 et 0,95 avec un pas de 0,2375.

Déformation des Matériaux

Ce cas d'usage porte 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é. On décrit 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 barre considérée est de longueur unitaire. Ce cas d'usage comporte une condition aux limites pour la contrainte de surface $t$, soit la quantité de travail nécessaire pour étirer la barre.

Les arguments des équations différentielles de MD sont définis sur une grille fixe, comme suit :

• $x$ est compris entre 0 et 1 avec un pas de 0,04.

Entrée

Pour exécuter la fonction Qiskit QUICK-PDE, tu peux ajuster les paramètres suivants :

Nom	Type	Description	Spécifique au cas d'usage	Exemple
use_case	Literal["MD", "CFD"]	Sélectionne le système d'équations différentielles à résoudre	Non	"CFD"
parameters	dict[str, Any]	Paramètres de l'équation différentielle (voir le tableau suivant pour plus de détails)	Non	{"a": 1.0, "b": 1.0}
nb_qubits	Optional[dict[str, dict[str, int]]]	Nombre de qubits par fonction et par variable. Des valeurs optimisées sont choisies par la fonction, mais si tu veux essayer de trouver une meilleure combinaison, tu peux remplacer les valeurs par défaut	Non	{"u": {"x": 1, "t":3}}
depth	Optional[dict[str, int]]	Profondeur de l'ansatz par fonction. Des valeurs optimisées sont choisies par la fonction, mais si tu veux essayer de trouver une meilleure combinaison, tu peux remplacer les valeurs par défaut	Non	{"u": 4}
optimizer	Optional[list[str]]	Optimiseurs à utiliser, soit "CMAES" de la bibliothèque Python cma, soit l'un des optimiseurs de scipy	MD	"SLSQP"
shots	Optional[list[int]]	Nombre de shots utilisés pour exécuter chaque circuit. Comme plusieurs étapes d'optimisation sont nécessaires, la longueur de la liste doit être égale au nombre d'optimiseurs utilisés (4 dans le cas CFD). Par défaut [50_000] * nb_optimizers pour MD et [5_00, 2_000, 5_000, 10_000] pour CFD	Non	[15_000, 30_000]
optimizer_options	Optional[dict[str, Any]]	Options à transmettre à l'optimiseur. Le détail de cette entrée dépend de l'optimiseur utilisé ; pour les spécificités, consulte la documentation de l'optimiseur concerné	MD	{"maxiter": 50 }
initialization	Optional[Literal["RANDOM", "PHYSICALLY_INFORMED"]]	Indique s'il faut démarrer avec des angles aléatoires ou des angles choisis intelligemment. Attention : les angles choisis intelligemment peuvent ne pas fonctionner dans 100 % des cas. Par défaut "PHYSICALLY_INFORMED".	Non	"RANDOM"
backend_name	Optional[str]	Le nom du backend à utiliser.	Non	"ibm_torino"
mode	Optional[Literal["job", "session", "batch"]]	Le mode d'exécution à utiliser. Par défaut "job".	Non	"job"

Les paramètres de l'équation différentielle (paramètres physiques et condition aux limites) doivent respecter le format suivant :

Cas d'usage	Clé	Type de valeur	Description	Exemple
CFD	a	float	Coefficient des valeurs initiales de $u$	1.0
CFD	b	float	Décalage des valeurs initiales de $u$	1.0
MD	t	float	contrainte de surface	12.0
MD	K	float	module de compressibilité	100.0
MD	n	int	loi de puissance	4.0
MD	b	float	force par unité de masse	10.0
MD	epsilon_0	float	limite de contrainte proportionnelle	0.1
MD	sigma_0	float	limite de déformation proportionnelle	5.0

Sortie

La sortie est un dictionnaire contenant la liste des points d'échantillonnage ainsi que les valeurs des fonctions en chacun de ces points :

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit-ibm-catalog

from numpy import array

solution = {
    "functions": {
        "u": array(
            [
                [0.01, 0.1, 1],
                [0.02, 0.2, 2],
                [0.03, 0.3, 3],
                [0.04, 0.4, 4],
            ]
        ),
    },
    "samples": {
        "t": array([0.1, 0.2, 0.3, 0.4]),
        "x": array([0.5, 0.6, 0.7]),
    },
}

La forme d'un tableau de solution dépend des échantillons des variables :

assert len(solution["functions"]["u"].shape) == len(solution["samples"])
for col_size, samples in zip(
    solution["functions"]["u"].shape, solution["samples"].values()
):
    assert col_size == len(samples)

La correspondance entre les points d'échantillonnage des variables de la fonction et les dimensions du tableau de solution est établie dans l'ordre alphanumérique du nom de la variable. Par exemple, si les variables sont "t" et "x", une ligne de solution["functions"]["u"] représente les valeurs de la solution pour un "t" fixé, et une colonne de solution["functions"]["u"] représente les valeurs de la solution pour un "x" fixé.

Voici un exemple illustrant comment obtenir la valeur de la fonction pour un ensemble de coordonnées spécifiques :

# 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

Benchmarks

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

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

Prise en main

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 :

from qiskit_ibm_catalog import QiskitFunctionsCatalog

catalog = QiskitFunctionsCatalog(channel="ibm_quantum_platform")

quick = catalog.load("colibritd/quick-pde")

Exemples

Pour commencer, essaie l'un des exemples suivants :

Dynamique des Fluides Numérique (CFD)

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", physical_parameters={"a": 1.0, "b": 0.0})

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

print(job.status())
solution = job.result()

'QUEUED'

import numpy as np
import matplotlib.pyplot as plt

_ = plt.figure()
ax = plt.axes(projection="3d")

# plot the solution using the 3d plotting capabilities of pyplot
t, x = np.meshgrid(solution["samples"]["t"], solution["samples"]["x"])
ax.plot_surface(
    t,
    x,
    solution["functions"]["u"],
    edgecolor="royalblue",
    lw=0.25,
    rstride=26,
    cstride=26,
    alpha=0.3,
)
ax.scatter(t, x, solution["functions"]["u"], marker=".")
ax.set(xlabel="t", ylabel="x", zlabel="u(t,x)")

plt.show()

Déformation des Matériaux

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

import matplotlib.pyplot as plt

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

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 faciliter le débogage, comme suit :

job = quick.run(use_case="mdf", physical_params={})

print(job.error_message())

# or write a wrapper around it for a more human readable version
def pprint_error(job):
    print("".join(eval(job.error_message())["error"]))

print("___")
pprint_error(job)

{"error": ["qiskit.exceptions.QiskitError: 'Unknown argument \"physical_params\", did you make a typo? -- https://docs.quantum.ibm.com/errors#1804'\n"]}
___
qiskit.exceptions.QiskitError: 'Unknown argument "physical_params", did you make a typo? -- https://docs.quantum.ibm.com/errors#1804'

Obtenir de l'aide

Pour toute assistance, contacte qiskit-function-support@colibritd.com.

Prochaines étapes

Recommandations

• Remplis le formulaire pour demander l'accès à la fonction QUICK-PDE.
• Essaie de modéliser l'écoulement d'un fluide non visqueux 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.