Skip to main content

View on GitHub

Open this notebook in GitHub to run it yourself
The Poisson equation is a partial differential equation that appears in various research fields, such as physics and engineering. It describes the distribution of a potential field, uu, under some source term bb: ∇2u=b,\nabla^2 u = b, where the ∇2\nabla^2 is the Laplacian (second derivatives) operator. One approach for numerically solving the Poisson equation is to move from the continuous description to a discrete one, using the finite difference method that casts the problem into a set of linear equations. Then, the solution can be obtained by a linear solver. In this notebook we treat the Poisson equation on a rectangular geometry, Lx×LyL_x\times L_y, with a Dirichlet boundary condition on the xx axis and a Neumann boundary condition on the yy axis: u(0)=u(Lx)=f0,    ∂yu∣y=0=∂yu∣y=Ly=g0.u(0) = u(L_x)=f_0,\,\,\,\, \partial_y u|_{y=0} = \partial_y u|_{y=L_y} = g_0. Furthermore, we assume that f0=g0=0f_0=g_0=0. The discretization of space, including the treatment of the above boundary conditions, is given in Figure
  1. The resulting linear equation reads:
L⋅u⃗=b⃗,       L=Lxx⊗Iy+Ix⊗Lyy,\mathcal{L}\cdot \vec{u} = \vec{b}, \,\,\,\,\,\,\, \mathcal{L} = \mathcal{L}_{xx} \otimes I_y + I_x \otimes \mathcal{L}_{yy}, where Lxx\mathcal{L}_{xx} and Lyy\mathcal{L}_{yy} are the Laplacian operators [1]: Lxx=1Δx2(3−10⋯0−12−1⋯00−12⋯0⋮⋮⋮⋱−100⋯−13),   Lyy=1Δy2(1−10⋯0−12−1⋯00−12⋯0⋮⋮⋮⋱−100⋯−11)\mathcal{L}_{xx} = \frac{1}{\Delta x^2} \begin{pmatrix} 3 & -1 & 0 & \cdots & 0 \\ -1 & 2 & -1 & \cdots & 0 \\ 0 & -1 & 2 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & -1 \\ 0 & 0 & \cdots & -1 & 3 \\ \end{pmatrix},\,\,\, \mathcal{L}_{yy} = \frac{1}{\Delta y^2} \begin{pmatrix} 1 & -1 & 0 & \cdots & 0 \\ -1 & 2 & -1 & \cdots & 0 \\ 0 & -1 & 2 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & -1 \\ 0 & 0 & \cdots & -1 & 1 \\ \end{pmatrix} and Δx\Delta x and Δy\Delta y are the discretization of the xx and yy axes, respectively. The square matrices above, which are of dimensions NxN_x and NyN_y, respectively, represent the solution at the inner part of our geometry.
In this notebook we solve the Poisson problem with the HHL quantum linear solver. We utilize similar ideas appearing in Ref.[2], where a quantum cosine and a quantum sine transforms [3] are performed towards achieving scalable implementation.

​
Building the Algorithm with Classiq

The HHL algorithm essentially applies a matrix inversion. Here we treat the Laplacian matrix, which can be diagonalized by quantum sine and cosine transforms. Thus, the matrix to invert is diagonal. The four main quantum blocks of the algorithm are thus (see Figure 2):
  1. Prepare the amplitudes of the source term on a quantum variable.
  2. Perform QST and QCT at the beginning of the computation.
This is done by applying the QST to the x qubits and the QCT to the y qubits.
  1. Perform matrix inversion for a diagonal matrix.
  2. Uncompute the QST and QCT at the end of the computation.
Screenshot 2025-07-01 at 0.10.30.png
Below we define several classical and quantum functions for constructing our quantum linear solver. 
import matplotlib.pyplot as plt
import numpy as np
import sympy

from classiq import *

​

​
Diagonal Hamiltonian Calculation

The eigenvalues of the Poisson equation with Dirichlet boundary conditions in the x direction and Neumann boundary conditions in the y direction are given by λk,j≡λx,k+λy,j\lambda_{k,j} \equiv \lambda_{x,k} +\lambda_{y,j} where λx,k=4Δx2sin⁡2(π2Nx(k+1)),        λy,j=4Δy2sin⁡2(π2Nyj),\lambda_{x,k} = \frac{4}{\Delta x^2} \sin^2\left(\frac{\pi}{2N_x} (k+1)\right),\,\,\,\,\,\,\,\, \lambda_{y,j} = \frac{4}{\Delta y^2} \sin^2\left(\frac{\pi}{2N_y} j\right), and k=0,1,…,Ny−1k = 0, 1, \ldots, N_y-1 and j=0,1,…,Nx−1j = 0, 1, \ldots, N_x-1. The HHL algorithm requires the application of an Hamiltonian simulation eiHe^{iH}. In this notebook this quantum block is implemented using the Suzuki-Trotter built-in function. We start with defining a function that gets NxN_x and NyN_y and returns the corresponding Hamiltonian. The decomposition of the diagonal matrix to the Pauli basis is done using the Walsh Hadamard transform. 
def get_poisson_dirichletx_neumanny_ham(nx, ny):
    dx2 = ny / nx
    dy2 = nx / ny
    eigenvalues_x = dx2 * 4 * np.sin(np.pi / 2 * np.arange(1, nx + 1) / nx) ** 2
    eigenvalues_y = dy2 * 4 * np.sin(np.pi / 2 * np.arange(0, ny) / ny) ** 2
    num_qubits_x = int(np.log2(nx))
    num_qubits_y = int(np.log2(ny))

    # Decompose the eigenvalues of the Poisson equation into Pauli terms and construct the corresponding Hamiltonian
    pauli_coefficients_x = sympy.fwht(eigenvalues_x / nx)
    pauli_coefficients_y = sympy.fwht(eigenvalues_y / ny)

    def convert_bitstring_to_pauli_representation(
        bitstring: int, num_qubits: int
    ) -> list[Pauli]:
        return [
            Pauli.Z if (bitstring >> s) & 1 else Pauli.I for s in range(num_qubits)
        ][::-1]

    hamiltonian_y = [
        PauliTerm(
            pauli=convert_bitstring_to_pauli_representation(i, num_qubits_y)
            + [Pauli.I] * num_qubits_x,
            coefficient=pauli_coefficients_y[i],
        )
        for i in range(ny)
    ]
    hamiltonian_x = [
        PauliTerm(
            pauli=[Pauli.I] * num_qubits_y
            + convert_bitstring_to_pauli_representation(i, num_qubits_x),
            coefficient=pauli_coefficients_x[i],
        )
        for i in range(nx)
    ]

    hamiltonian = hamiltonian_x + hamiltonian_y
    return hamiltonian

​

​
Hamiltonian Evolution for QPE

The HHL is based on a QPE applied on eiHte^{iHt}. For this, we need to define a function implementing (eiHt)p\left(e^{iHt}\right)^p for an integer power pp. Since in our case the Hamiltonian is diagonal, an exact implementation is given by the first order Suzuki-Trotter formula, where in addition (eiHt)p=eipHt\left(e^{iHt}\right)^p = e^{ipHt}. 
@qfunc
def powered_hamiltonian_evolution(
    hamiltonian: CArray[PauliTerm],  # the hamiltonian H
    scaling: CReal,  # the scaling factor t
    p: CInt,  # the power
    qba: QArray,
):
    suzuki_trotter(
        pauli_operator=hamiltonian,
        evolution_coefficient=p * (-2 * np.pi * scaling),
        order=1,
        repetitions=1,
        qbv=qba,
    )

​

​
Sine and Cosine Transforms

For the system treated in this notebook, diagonalization of the 2D2D Laplacian is given by appying a quantum sine (cosine) transform on the xx (yy) dimension. We define a quantum function for implementing this, using the open library functions: 
@qfunc
def qsct_2d(xy_variable: QArray[QNum, 2]):
    qst_type2(xy_variable[0])
    qct_type2(xy_variable[1])

​

​
Matrix Inversion

Finally, we define a matrix inversion using a QPE routine, similar to the example given in the basic HHL notebook. 
@qfunc
def matrix_inversion_HHL(
    prefactor: float,
    my_unitary: QCallable[CInt, QArray],
    state: QArray,
    phase: QNum,
    indicator: Output[QBit],
):
    allocate(indicator)

    within_apply(
        within=lambda: qpe_flexible(
            unitary_with_power=lambda power: my_unitary(power, state),
            phase=phase,
        ),
        apply=lambda: assign_amplitude_table(
            lookup_table(lambda p: 0 if p == 0 else prefactor / p, phase),
            phase,
            indicator,
        ),
    )

​
Example: Non-separable Source Term

We solve an example with a square grid of Nx, Ny=23N_x,\, N_y=2^3. For the source term we take a non-separable 2Nx+Ny2^{N_x+N_y} vector that represents the function b=F[xy(x−Lx)(y−Ly)].b = F\left[xy(x-L_x)(y-L_y)\right]. We choose F(z)=tanh⁡(z)2F(z)=\tanh(z)^2, which satisfies the boundary conditions. 
# Set discretization of X and Y axes.
NUM_QUBITS_X = 3
NUM_QUBITS_Y = 3
nx = 2**NUM_QUBITS_X
ny = 2**NUM_QUBITS_Y

# Set the source term as function of x and y such that b(x,y) = exp(-1/xy(x-Lx)(y-Ly))
xgrid = (np.arange(nx) + 0.5) / nx
ygrid = (np.arange(ny) + 0.5) / ny
zgrid = np.kron(ygrid, xgrid) * np.kron(
    (ygrid - 1), (xgrid - 1)
)  # in classiq the variable order is [y,x]
b_vector = np.tanh(zgrid) ** 2

# Normalize the source term
b_vector = b_vector / np.linalg.norm(b_vector)
We can plot the source term: 
b_matrix = b_vector.reshape((nx, ny))
xmesh, ymesh = np.meshgrid(xgrid, ygrid)
plt.contourf(xmesh, ymesh, b_matrix.transpose())
plt.xlabel("x")
plt.ylabel("y")
ax = plt.gca()
ax.axis("equal")
plt.colorbar()
plt.show()
output Next, we calculate the Hamiltonian for the specific discretization: 
hamiltonian = get_poisson_dirichletx_neumanny_ham(nx, ny)
Now, we need to define the resolution (QPE size) for the quantum solver. The required number of QPE phase qubits depends on the condition number κ\kappa of the inverted matrix. In the case of the Laplacian matrix, this parameter is of the order of the matrix dimension O(2Nx+Ny)O(2^{N_x+N_y}), which eliminates the exponential advantage. However, in our example we have a smooth source term, and high modes are expected to have minor effects on the solution. The highest mode of L\mathcal{L} is λmax⁡=8\lambda_{\max}=8, whereas the smallest one is λx,0∼4π2/(2Nx)2\lambda_{x,0}\sim 4\pi^2/(2N_x)^2. Let us assume that the highest mode participating in the solution is ∼26λx,0\sim 2^6 \lambda_{x,0}, and take a QPE size of six qubits. 
# number of qubits for the QPE
QPE_SIZE = 6
MAX_EIG = 4 * (np.pi / (2 * nx)) ** 2 * 2**6
# parameters for the amplitude preparation
PREFACTOR = 2**-QPE_SIZE
We build the model and synthesize it: Comment: the model is designed for the case of symmetric grid Nx=NyN_x=N_y, working with QArray[QNum,2]. In the more general case one can define a QStruct with two QNum variable of different sizes 
@qfunc
def main(
    xy_variable: Output[QArray[QNum[NUM_QUBITS_X], 2]],
    phase_var: Output[QNum[QPE_SIZE]],
    indicator: Output[QBit],
):

    prepare_amplitudes(b_vector.tolist(), 0.0, xy_variable)
    allocate(phase_var)

    within_apply(
        within=lambda: qsct_2d(xy_variable),
        apply=lambda: matrix_inversion_HHL(
            prefactor=PREFACTOR,
            my_unitary=lambda p, target: powered_hamiltonian_evolution(
                hamiltonian=hamiltonian,
                scaling=1 / MAX_EIG,
                p=p,
                qba=target,
            ),
            state=xy_variable,
            phase=phase_var,
            indicator=indicator,
        ),
    )


qmod = create_model(main, constraints=Constraints(max_width=18))
qprog = synthesize(qmod)
show(qprog)
Output:

Quantum program link: https://platform.classiq.io/circuit/32pZbY4ZtlxYiymOopvtvjloFFz

​

​
Executing and Plotting the Result

We run the quantum program on a statevector simulator to retrieve the full solution. 
from classiq.execution import (
    ClassiqBackendPreferences,
    ClassiqSimulatorBackendNames,
    ExecutionDetails,
    ExecutionPreferences,
    set_quantum_program_execution_preferences,
)

execution_preferences = ExecutionPreferences(
    backend_preferences=ClassiqBackendPreferences(backend_name="simulator_statevector")
)
qprog = set_quantum_program_execution_preferences(qprog, execution_preferences)
result = execute(qprog).result_value()


import pandas as pd


# Retrieve the solution of the Poisson equation from the results of the quantum program.
def extract_result_from_statevector_simulation(
    df: pd.DataFrame,
    normalization: float,
    indicator_name: str,
    phase_name: str,
    xy_name: str,
) -> np.ndarray:

    # Initialize a matrix
    state_matrix = np.zeros((nx, ny), dtype=complex)

    # Filter only the successful states.
    filtered_st = df[
        (df[indicator_name] == 1)
        & (df[phase_name] == 0)
        & (np.abs(df.amplitude) > 1e-12)
    ]

    # Allocate values
    state_matrix[tuple(np.array(filtered_st[xy_name].tolist()).T)] = (
        filtered_st.amplitude.to_numpy()
    )

    # Normalize the solution of the Poisson equation.
    global_phase = np.angle(state_matrix[0, 0])
    result_matrix = np.real(state_matrix / np.exp(1j * global_phase)) / normalization

    return result_matrix


result_matrix = extract_result_from_statevector_simulation(
    result.dataframe, PREFACTOR / MAX_EIG, "indicator", "phase_var", "xy_variable"
)
The resulting statevector is given up to a sign. The sign in the center of the solution is expected to be positive. We can correct accordingly: 
if result_matrix[nx // 2, ny // 2] < 0:
    result_matrix = -result_matrix
Finally, we print the result and compare it to the classical solution: 

## getting the classical solution
Hx = (
    2 * np.diag(np.ones(nx))
    - np.diag(np.ones(nx - 1), 1)
    - np.diag(np.ones(nx - 1), -1)
)
Hx[0, 0] = Hx[nx - 1, nx - 1] = 3

Hy = (
    2 * np.diag(np.ones(ny))
    - np.diag(np.ones(ny - 1), 1)
    - np.diag(np.ones(ny - 1), -1)
)
Hy[0, 0] = Hy[ny - 1, ny - 1] = 1

# we need to flip the order of x and y
b_classical = (b_matrix.T).reshape(nx * ny)

H = np.kron(Hx, np.identity(ny)) + np.kron(np.identity(nx), Hy)
classical_result_vector = np.linalg.solve(H, b_classical)
classical_result_matrix = classical_result_vector.reshape((nx, ny))


xmesh, ymesh = np.meshgrid(xgrid, ygrid)

fig, axs = plt.subplots(1, 2, figsize=(12, 6))
contour0 = axs[0].contourf(xmesh, ymesh, result_matrix.transpose())
axs[0].axis("equal")
axs[0].axis("square")
axs[0].title.set_text("Quantum")
fig.colorbar(contour0, ax=axs[0])

contour1 = axs[1].contourf(xmesh, ymesh, classical_result_matrix.transpose())
axs[1].axis("equal")
axs[1].axis("square")
axs[1].title.set_text("Classical")
fig.colorbar(contour1, ax=axs[1])

plt.show()
output

​
References

[1]: Strang, G., 1999 SIAM Review 41 135. The discrete cosine transform. [2]: Yudong Cao et al., 2013 New J. Phys. 15 013021. Quantum algorithm and circuit design solving the Poisson equation. [3]: Klappenecker, A., & Rotteler M., Discrete cosine transforms on quantum computers.