QUICK-PDE: Una Qiskit Function de ColibriTD

Consulta la referencia de la API

nota

Las Qiskit Functions son una funcionalidad experimental disponible para los usuarios de los planes IBM Quantum® Premium, Flex y On-Prem (mediante la API de IBM Quantum Platform). Se encuentran en estado de versión preliminar y están sujetas a cambios.

Descripción general

El solucionador de Ecuaciones Diferenciales Parciales (PDE) que se presenta aquí forma parte de nuestra plataforma Quantum Innovative Computing Kit (QUICK) (QUICK-PDE) y está empaquetado como una Qiskit Function. Con la función QUICK-PDE puedes resolver ecuaciones diferenciales parciales específicas de dominio en QPUs de IBM Quantum. Esta función se basa en el algoritmo descrito en el artículo de descripción del H-DES de ColibriTD. Este algoritmo puede resolver problemas complejos de múltiples físicas, comenzando con Dinámica de Fluidos Computacional (CFD) y Deformación de Materiales (MD), con más casos de uso próximamente.

Para abordar las ecuaciones diferenciales, las soluciones de prueba se codifican como combinaciones lineales de funciones ortogonales (típicamente polinomios de Chebyshev, y más específicamente $2^n$ de ellos, donde $n$ es el número de qubits que codifican tu función), parametrizadas por los ángulos de un Circuito Cuántico Variable (VQC). El ansatz genera un estado que codifica la función, que se evalúa mediante observables cuyas combinaciones permiten evaluar la función en todos los puntos. Luego puedes evaluar la función de pérdida en la que se codifican las ecuaciones diferenciales y ajustar los ángulos en un bucle híbrido, como se muestra a continuación. Las soluciones de prueba se acercan gradualmente a las soluciones reales hasta alcanzar un resultado satisfactorio.

Además de este bucle híbrido, también puedes encadenar distintos optimizadores. Esto es útil cuando quieres que un optimizador global encuentre un buen conjunto de ángulos y luego un optimizador más fino siga un gradiente hasta el mejor conjunto de ángulos vecinos. En el caso de la dinámica de fluidos computacional (CFD), la secuencia de optimización predeterminada produce los mejores resultados; en el caso de la deformación de materiales (MD), aunque la opción predeterminada proporciona buenos resultados, puedes configurarla con mayor detalle para obtener beneficios específicos del problema.

Ten en cuenta que para cada variable de la función especificamos el número de qubits (con el que puedes experimentar). Al apilar 10 circuitos idénticos y evaluar los 10 observables idénticos en diferentes qubits a lo largo de un gran circuito, puedes aplicar mitigación de ruido dentro del proceso de optimización CMA, apoyándote en el método del noise learner, y reducir significativamente el número de shots necesarios.

Dinámica de fluidos computacional

La ecuación de Burgers no viscosa modela fluidos no viscosos en movimiento de la siguiente manera:

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

$u$ representa el campo de velocidad del fluido. Este caso de uso tiene una condición de frontera temporal: puedes seleccionar la condición inicial y luego permitir que el sistema se relaje. Actualmente, las únicas condiciones iniciales aceptadas son funciones lineales: $ax + b$. La solución analítica es:

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

Las ecuaciones de Euler sin presión modelan el flujo de fluido no viscoso compresible con amortiguamiento de la siguiente manera:

$\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$ representa el campo de densidad, $u$ el campo de velocidad y $\mu$ un coeficiente de amortiguamiento. En nuestra formulación, establecemos $\lambda = 1$, por lo que no se usará como parámetro en lo que sigue. Este caso de uso tiene condiciones de contorno temporales: $g(0, x) = e^{-x}$ y $u(0, x) = x$. La solución analítica es:

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

Los argumentos para las ecuaciones diferenciales de CFD se encuentran en una cuadrícula fija, de la siguiente manera:

• $t$ está entre 0 y 0.95 con 41 puntos de muestreo. $x$ está entre 0 y 0.95 con 41 puntos de muestreo.

Deformación de materiales

Este caso de uso se centra en la deformación hipoelástica con el ensayo de tracción unidimensional, en el que una barra fija en el espacio es jalada por su otro extremo. Describimos el problema de la siguiente manera:

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

$\sigma' - b = 0,$

$K$ representa el módulo volumétrico del material que se estira, $n$ el exponente de una ley de potencias, $b$ la fuerza por unidad de masa, $\epsilon_0$ el límite de tensión proporcional, $\sigma_0$ el límite de deformación proporcional, $u$ la función de tensión y $\sigma$ la función de deformación. La solución analítica es:

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

donde $\sigma_0 = g(0)$ es la condición de contorno de la deformación en $x=0$.

La barra considerada tiene longitud unitaria. Este caso de uso tiene una condición de frontera para la tensión superficial $t$, es decir, la cantidad de trabajo necesaria para estirar la barra.

Los argumentos para las ecuaciones diferenciales de MD se encuentran en una cuadrícula fija, de la siguiente manera:

• $x$ está entre 0 y 1 con 30 puntos de muestreo.

Benchmarks

La siguiente tabla presenta estadísticas de varias ejecuciones de nuestra función.

Ejemplo	Número de qubits	Inicialización	Error	Tiempo total (min)	Uso del runtime (min)
Ecuación de Burgers no viscosa	50	PHYSICALLY_INFORMED	$10^{-2}$	66	25
Ecuaciones de Euler sin presión	70	PHYSICALLY_INFORMED	$10^{-2}$	48	34
Ensayo de tracción 1D hipoelástico	18	RANDOM	$10^{-2}$	123	100

Comenzar

Rellena el formulario para solicitar acceso a la función QUICK-PDE. Luego, suponiendo que ya has guardado tu cuenta en tu entorno local, selecciona la función de la siguiente manera:

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

Ejemplos

Para empezar, prueba uno de los siguientes ejemplos:

Ecuación de Burgers no viscosa (CFD)

Para la ecuación de Burgers, cuando las condiciones iniciales se establecen como $u(0,x) = x$, los resultados son los siguientes:

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

Verifica el estado de la carga de trabajo de tu Qiskit Function u obtén los resultados de la siguiente manera:

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

Ecuación de Euler sin presión (CFD)

Para la ecuación de Euler, cuando las condiciones iniciales se establecen como $g(0, x) = e^{-x}$ y $u(0, x) = x$, para un $\mu$ dado (aquí, $\mu = 0.1$) y $\lambda = 1$, los resultados son los siguientes:

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

Deformación de materiales

El caso de uso de deformación de materiales requiere los parámetros físicos del material y la fuerza aplicada, de la siguiente manera:

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

El siguiente es un ejemplo de cómo obtener el valor de la función para un conjunto específico de coordenadas:

# 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

Obtener mensajes de error

Si el estado de tu carga de trabajo es ERROR, usa job.error_message() para obtener el mensaje de error que te ayude a depurar, de la siguiente manera:

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)

Obtener soporte

Para obtener soporte, contacta a qiskit-function-support@colibritd.com.

Próximos pasos

Recomendaciones

• Rellena el formulario para solicitar acceso a la función QUICK-PDE.
• Visita la referencia de la API para esta Qiskit Function.
• Prueba modelar un fluido no viscoso en movimiento con QUICK-PDE en el tutorial.
• Consulta Jaffali, H., et al. (2025). H-DES: a Quantum-Classical Hybrid Differential Equation Solver. arXiv preprint arXiv:2410.01130.