# Burton-Miller formulation is not correct

I found that the Burton-Miller formulation yields higher errors compared to the standard boundary element equation in my test case. Can someone explain why this might be happening?

Reproduction code:

``````import bempp.api
import numpy as np

bempp.api.enable_console_logging("debug")
bempp.api.BOUNDARY_OPERATOR_DEVICE_TYPE = "gpu"
bempp.api.POTENTIAL_OPERATOR_DEVICE_TYPE = "gpu"

x0 = np.array([0.0, 0.0, 0.0])
rerr_lst = []
wave_number = 3

@bempp.api.complex_callable
def get_dirichlet(x, n, domain_index, result):
r_vec = x - x0
r = np.dot(r_vec, r_vec) ** 0.5
k = wave_number
result[0] = np.exp(1j * k * r) / (4 * np.pi * r)

@bempp.api.complex_callable
def get_neumann(x, n, domain_index, result):
r_vec = x - x0
r = np.dot(r_vec, r_vec) ** 0.5
k = wave_number
ikr = 1j * k * r
result[0] = (
-np.exp(ikr) / (4 * np.pi * r * r * r) * (1 - ikr) * ((x - x0) * n).sum()
)

wavelength = 2 * np.pi / abs(wave_number)
h = wavelength / 6
grid = bempp.api.shapes.ellipsoid(3, 1, 1, h=h * 0.5)

space = bempp.api.function_space(grid, "P", 1)

identity = bempp.api.operators.boundary.sparse.identity(space, space, space)
slp = bempp.api.operators.boundary.helmholtz.single_layer(
space, space, space, wave_number
)
dlp = bempp.api.operators.boundary.helmholtz.double_layer(
space, space, space, wave_number
)
hyp = bempp.api.operators.boundary.helmholtz.hypersingular(
space, space, space, wave_number
)
space, space, space, wave_number
)

LHS = 0.5 * identity - dlp
RHS = -slp

neumann_fun = bempp.api.GridFunction(space, fun=get_neumann)

RHS = RHS * neumann_fun

dirichlet_fun, info = bempp.api.linalg.gmres(LHS, RHS, tol=1e-5)
dirichlet_fun_gt = bempp.api.GridFunction(space, fun=get_dirichlet)

rerr = np.linalg.norm(
dirichlet_fun.coefficients - dirichlet_fun_gt.coefficients
) / np.linalg.norm(dirichlet_fun_gt.coefficients)

rerr_lst.append(rerr)

beta = 1j / wave_number
LHS = 0.5 * identity - dlp - beta * hyp
RHS = (-slp - beta * (adlp + 0.5 * identity)) * neumann_fun

dirichlet_fun, info = bempp.api.linalg.gmres(LHS, RHS, tol=1e-5)
dirichlet_fun_gt = bempp.api.GridFunction(space, fun=get_dirichlet)

rerr = np.linalg.norm(
dirichlet_fun.coefficients - dirichlet_fun_gt.coefficients
) / np.linalg.norm(dirichlet_fun_gt.coefficients)

rerr_lst.append(rerr)
print(rerr_lst)

``````

The output is

``````[0.00011713474983450146, 1.5905369315737763]
``````

The error is too high, I even wonder if the Bempp implementation is wrong?

I found the problem. There is a negative sign in front of the hypersingular operator defined in Bempp, so negating the hypersingular operator solves the issue.