QUICK-PDE: Uma Função Qiskit da ColibriTD

Veja a referência de API

nota

As Qiskit Functions são um recurso experimental disponível para usuários dos planos IBM Quantum® Premium, Flex e On-Prem (via API da IBM Quantum Platform). Elas estão em status de versão prévia e sujeitas a alterações.

Visão geral

O solver de Equações Diferenciais Parciais (EDP) apresentado aqui faz parte da nossa plataforma Quantum Innovative Computing Kit (QUICK) (QUICK-PDE), e é disponibilizado como uma Função Qiskit. Com a função QUICK-PDE, você pode resolver equações diferenciais parciais específicas de domínio em QPUs da IBM Quantum. Esta função é baseada no algoritmo descrito no artigo de descrição do H-DES da ColibriTD. Esse algoritmo pode resolver problemas multifísicos complexos, começando com Dinâmica de Fluidos Computacional (DFC) e Deformação de Materiais (DM), com outros casos de uso chegando em breve.

Para lidar com as equações diferenciais, as soluções tentativas são codificadas como combinações lineares de funções ortogonais (tipicamente polinômios de Chebyshev, e mais especificamente $2^n$ deles, onde $n$ é o número de qubits que codificam sua função), parametrizadas pelos ângulos de um Circuito Quântico Variacional (VQC). O ansatz gera um estado que codifica a função, avaliado por observáveis cujas combinações permitem calcular a função em todos os pontos. Você pode então avaliar a função de perda na qual as equações diferenciais são codificadas e ajustar os ângulos em um laço híbrido, como mostrado a seguir. As soluções tentativas se aproximam gradualmente das soluções reais até que você alcance um resultado satisfatório.

Além desse laço híbrido, você também pode encadear diferentes otimizadores. Isso é útil quando você quer que um otimizador global encontre um bom conjunto de ângulos e, em seguida, um otimizador mais refinado siga um gradiente até o melhor conjunto de ângulos vizinhos. No caso de dinâmica de fluidos computacional (DFC), a sequência de otimização padrão produz os melhores resultados — mas no caso de deformação de materiais (DM), embora o padrão forneça bons resultados, você pode ajustá-lo ainda mais para obter ganhos específicos do problema.

Observe que, para cada variável da função, especificamos o número de qubits (com o qual você pode experimentar). Ao empilhar 10 circuitos idênticos e avaliar os 10 observáveis idênticos em qubits diferentes ao longo de um único circuito grande, você pode mitigar ruído dentro do processo de otimização CMA, com base no método de aprendizado de ruído, e reduzir significativamente o número de shots necessários.

Dinâmica de fluidos computacional

A equação de Burgers invíscida modela fluidos não viscosos em escoamento da seguinte forma:

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

$u$ representa o campo de velocidade do fluido. Esse caso de uso tem uma condição de contorno temporal: você pode selecionar a condição inicial e então permitir que o sistema relaxe. Atualmente, as únicas condições iniciais aceitas são funções lineares: $ax + b$. A solução analítica é:

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

As equações de Euler sem pressão modelam o fluxo de fluido invíscido compressível com amortecimento da seguinte forma:

$\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 o campo de densidade, $u$ o campo de velocidade e $\mu$ um coeficiente de amortecimento. Em nossa formulação, definimos $\lambda = 1$, portanto ele não será usado como parâmetro a seguir. Esse caso de uso tem condições de contorno temporais: $g(0, x) = e^{-x}$ e $u(0, x) = x$. A solução analítica é:

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

Os argumentos para as equações diferenciais da DFC estão em uma grade fixa, como segue:

• $t$ está entre 0 e 0,95 com 41 pontos de amostra. $x$ está entre 0 e 0,95 com 41 pontos de amostra.

Deformação de Materiais

Este caso de uso foca na deformação hipoelástica com o teste de tração unidimensional, no qual uma barra fixada no espaço é puxada em sua outra extremidade. Descrevemos o problema da seguinte forma:

$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 o módulo volumétrico do material sendo esticado, $n$ o expoente de uma lei de potência, $b$ a força por unidade de massa, $\epsilon_0$ o limite de tensão proporcional, $\sigma_0$ o limite de deformação proporcional, $u$ a função de tensão e $\sigma$ a função de deformação. A solução analítica é:

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

onde $\sigma_0 = g(0)$ é a condição de contorno na deformação em $x=0$.

A barra considerada tem comprimento unitário. Este caso de uso tem uma condição de contorno para a tração superficial $t$, ou seja, a quantidade de trabalho necessária para esticar a barra.

Os argumentos para as equações diferenciais da DM estão em uma grade fixa, como segue:

• $x$ está entre 0 e 1 com 30 pontos de amostra.

Benchmarks

A tabela a seguir apresenta estatísticas sobre várias execuções da nossa função.

Example	Number of qubits	Initialization	Error	Total time (min)	Runtime usage (min)
Inviscid Burgers' equation	50	PHYSICALLY_INFORMED	$10^{-2}$	66	25
Pressureless Euler's equations	70	PHYSICALLY_INFORMED	$10^{-2}$	48	34
Hypoelastic 1D tensile test	18	RANDOM	$10^{-2}$	123	100

Primeiros passos

Preencha o formulário para solicitar acesso à função QUICK-PDE. Em seguida, assumindo que você já salvou sua conta em seu ambiente local, selecione a função da seguinte forma:

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

Exemplos

Para começar, experimente um dos seguintes exemplos:

Equação de Burgers invíscida (DFC)

Para a equação de Burgers, quando as condições iniciais são definidas como $u(0,x) = x$, os resultados são os seguintes:

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

Verifique o status ou retorne os resultados da carga de trabalho da sua Função Qiskit da seguinte forma:

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

Equação de Euler sem pressão (CFD)

Para a equação de Euler, quando as condições iniciais são definidas como $g(0, x) = e^{-x}$ e $u(0, x) = x$, para um dado $\mu$ (aqui, $\mu = 0.1$) e $\lambda = 1$, os resultados são os seguintes:

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

Deformação de Materiais

O caso de uso de deformação de materiais requer os parâmetros físicos do seu material e a força aplicada, da seguinte forma:

# 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 seguir há um exemplo de como obter o valor da função para um 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

Recuperar mensagens de erro

Se o status da sua carga de trabalho for ERROR, use job.error_message() para recuperar a mensagem de erro e ajudar na depuração, da seguinte forma:

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)

Obter suporte

Para suporte, entre em contato com qiskit-function-support@colibritd.com.

Próximos passos

Recomendações

• Preencha o formulário para solicitar acesso à função QUICK-PDE.
• Visite a referência de API para esta Função Qiskit.
• Experimente modelar um fluido não viscoso em escoamento usando QUICK-PDE no tutorial.
• Consulte Jaffali, H., et al. (2025). H-DES: a Quantum-Classical Hybrid Differential Equation Solver. arXiv preprint arXiv:2410.01130.