QUICK-PDE: Una función Qiskit de ColibriTD

Consulte la referencia de la API

Visión general

El solucionador de ecuaciones diferenciales parciales (PDE) que aquí se presenta forma parte de nuestra plataforma Quantum Innovative Computing Kit (QUICK) (QUICK-PDE) y se incluye como una función Qiskit. Con la función QUICK-PDE, puede resolver ecuaciones diferenciales parciales específicas de dominio en las QPU de IBM Quantum. Esta función se basa en el algoritmo descrito en el documento descriptivo de H-DES de ColibriTD's. Este algoritmo puede resolver problemas multifísicos complejos, comenzando por la dinámica de fluidos computacional (CFD) y la deformación de materiales (MD), y pronto se añadirán otros casos de uso.

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 concretamente $2^n$ de ellos donde $n$ es el número de qubits que codifican su 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. A continuación, puede 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 que se alcanza un resultado satisfactorio.

Además de este bucle híbrido, también se pueden encadenar distintos optimizadores. Esto es útil cuando se desea un optimizador global para encontrar un buen conjunto de ángulos, y luego un optimizador más afinado para seguir un gradiente al 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, pero en el caso de la deformación de materiales (MD), aunque la predeterminada proporciona buenos resultados, puede configurarla aún más para obtener ventajas específicas del problema.

Observa que para cada variable de la función, especificamos el número de qubits (con el que puedes jugar). Al apilar 10 circuitos idénticos y evaluar los 10 observables idénticos en diferentes qubits a lo largo de un gran circuito, se puede mitigar el ruido dentro del proceso de optimización CMA, confiando en el método del aprendiz de ruido, y reducir significativamente el número de disparos necesarios.

Dinámica de fluidos computacional

La ecuación de Burgers para fluidos no viscosos modela el flujo de dichos fluidos de la siguiente manera:

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

$u$ representa el campo de velocidades del fluido. Este caso de uso tiene una condición límite temporal: puedes seleccionar la condición inicial y, a continuación, dejar que el sistema se estabilice. 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 un fluido compresible e inviscido 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 amortiguación. En nuestra formulación, establecemos que $\lambda = 1$, por lo que no se utilizará 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 de las ecuaciones diferenciales de CFD están en una rejilla fija, como sigue:

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

Deformación del material

Este caso de uso se centra en la deformación hipoelástica mediante un ensayo de tracción unidimensional, en el que se aplica una fuerza de tracción en un extremo de una barra fijada en el espacio. 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 de compresibilidad del material sometido a estiramiento, $n$ el exponente de una ley de potencia, $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 relativa a la deformación en $x=0$.

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

Los argumentos de las ecuaciones diferenciales de MD están en una rejilla fija, como sigue:

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

Referencias comparativas

El siguiente cuadro presenta estadísticas sobre varias ejecuciones de nuestra función.

Cómo empezar

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

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 Inviscid Burgers (CFD)

En el caso de la ecuación de Burgers, cuando las condiciones iniciales se establecen en $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}
)

Comprueba el estado de tu carga de trabajo de Qiskit Function o 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)

En el caso de la ecuación de Euler, cuando las condiciones iniciales se establecen en $g(0, x) = e^{-x}$ y $u(0, x) = x$, para un $\mu$ dado (en este caso, $\mu = 0.1$ ) y un $\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 del material

El caso de uso de la deformación del material requiere los parámetros físicos de su material y la fuerza aplicada, como se indica a continuación:

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

A continuación se muestra un ejemplo de cómo obtener el valor de la función para un conjunto concreto 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 su carga de trabajo es ERROR, utilice job.error_message() para obtener el mensaje de error para ayudar 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 asistencia, póngase en contacto con [email protected].

Próximos pasos