# Boundary element method (BEM) in Mathematica

Please is there any mathematica codes for solving PDEs say Laplace or Poisson Equations by using boundary elements method?

Best regards,

• There is no build-in support for BEM in Mathematica as far as I know. There is this old package library.wolfram.com/infocenter/Articles/3802 we introduce two new Mathematica packages for solving problems arising from the engineering field, by using the Boundary Elements Method (BEM) Mar 17, 2020 at 20:46
• Dear Nasser, many thanks for your super fast response, how can I got the codes for above? Mar 17, 2020 at 20:54
• Sorry, I do not know. it is old as you can see, from 1997. May be you contact the author at personales.unican.es/iglesias or wait to see if someone knows of more recent packages. It used to be that library.worlfram will have the software or paper there on same page, but I do not think this site is being maintained for long time now. Mar 17, 2020 at 21:12
• The article is available at witpress.com/elibrary/wit-transactions-on-engineering-sciences/… Mar 17, 2020 at 23:02
• This would be a nice addition; if you get around to write / find a package I'd very much like to see it. I'd appreciate if you could @ ping me. Mar 18, 2020 at 7:05

I have a working code that I translated from MATLAB code. I improved it a little, but I didn’t finish it to the end. I give it as it is. Code BEM to solve the Laplace equation in a rectangle

N0 = NO = 20;(*the number of boundary elements per side*)
NN = 4*NO;
dl = 1/NO;
Do[ xb[i] = (i - 1)*dl;
yb[i] = 0;
xb[NO + i] = 1;
yb[NO + i] = xb[i];
xb[2*NO + i] = 1 - xb[i];
yb[2 NO + i] = 1;
xb[3*NO + i] = 0;
yb[3*NO + i] = 1 - xb[i];, {i, 1, NO}];

xb[NN + 1] = xb[1];
yb[NN + 1] = yb[1];

Do[xm[i] = 0.5*(xb[i] + xb[i + 1]);
ym[i] = 0.5*(yb[i] + yb[i + 1]);
lg[i] = Sqrt[(xb[i + 1] - xb[i])^2 + (yb[i + 1] - yb[i])^2];
nx[i] = (yb[i + 1] - yb[i])/lg[i];
ny[i] = (-xb[i + 1] + xb[i])/lg[i];, {i, 1, NN}];

(*boundary conditions*)
Do[BCT[i] = 0;
BCV[i] = 0;, {i, 1, NN}];

Do[BCT[i] = 1;
BCV[i] = 0;, {i, 1, N0}];

Do[BCT[i] = 0;
BCV[i] = Cos[Pi*ym[i]];, {i, N0 + 1, 2 N0}];

Do[BCT[i] = 1;
BCV[i] = 0;, {i, 2 N0 + 1, 3 N0}];

(*end*)
Do[BC[m] = 0;
Do[A[k] = lg[k]^2;
B[k] = 2*lg[k]*(-ny[k]*(xb[k] - xm[m]) + nx[k]*(yb[k] - ym[m]));
E0[k] = (xb[k] - xm[m])^2 + (yb[k] - ym[m])^2;
D0[k] = Sqrt[Abs[4*A[k]*E0[k] - B[k]^2]];
BA[k] = B[k]/A[k];
EA[k] = E0[k]/A[k];
If[ 4*A[k]*E0[k] - B[k]^2 == 0,
PF1[k] =
0.5*lg[k]*(Log[lg[k]] + (1 + 0.5*BA[k])*Log[Abs[1 + 0.5*BA[k]]] -
0.5*BA[k]*Log[Abs[0.5*BA[k]]] - 1),
PF1[k] =
0.25*lg[k]*(2*(Log[lg[k]] - 1) -
0.5*BA[k]*Log[Abs[EA[k]]] + (1 + 0.5*BA[k])*
Log[Abs[1 + BA[k] + EA[k]]] + (D0[k]/
A[k])*(ArcTan[D0[k], 2*A[k] + B[k]] - ArcTan[D0[k], B[k]]))];
If[ 4*A[k]*E0[k] - B[k]^2 == 0, PF2[k] = 0,
PF2[k] =
lg[k]*(nx[k]*(xb[k] - xm[m]) + ny[k]*(yb[k] - ym[m]))/
D0[k]*(ArcTan[D0[k], 2*A[k] + B[k]] - ArcTan[D0[k], B[k]])];
F1[k] = PF1[k]/Pi;
F2[k] = PF2[k]/Pi;

If[ k == m,
del = 1, del = 0];
If[ BCT[k] == 0,
AB[m, k] = -F1[k], AB[m, k] = F2[k] - 0.5*del];
If[ BCT[k] == 0, BC[m] = BC[m] + BCV[k]*(-F2[k] + 0.5*del),
BC[m] = BC[m] + BCV[k]*F1[k]];, {k, 1, NN}];, {m, 1, NN}]

NAC = NN;

ABM = SparseArray[{i_, j_} -> AB[i, j], {NN, NN}];

BCM = SparseArray[{i_} -> BC[i], {NN}]; Z = LinearSolve[ABM, BCM];

Do[
If [BCT[m] == 0,
phi[m] = BCV[m], phi[m] = Z[[m]]];
If[ BCT[m] == 0, dphi[m] = Z[[m]],
dphi[m] = BCV[m]];, {m, 1, NN}];

sol = ParallelTable[sum = 0; Do[A[i] = lg[i]^2;
B[i] = 2*lg[i]*(-ny[i]*(xb[i] - xi) + nx[i]*(yb[i] - eta));
E1[i] = (xb[i] - xi)^2 + (yb[i] - eta)^2;
D1[i] = Sqrt[Abs[4*A[i]*E1[i] - B[i]^2]];
BA[i] = B[i]/A[i];
EA[i] = E1[i]/A[i];
If[ 4*A[i]*E1[i] - B[i]^2 == 0,
PF1[i] =
0.5*lg[i]*(Log[lg[i]] + (1 + 0.5*BA[i])*
Log[Abs[1 + 0.5*BA[i]]] - 0.5*BA[i]*Log[Abs[0.5*BA[i]]] -
1), PF1[i] =
0.25*lg[i]*(2*(Log[lg[i]] - 1) -
0.5*BA[i]*Log[Abs[EA[i]]] + (1 + 0.5*BA[i])*
Log[Abs[1 + BA[i] + EA[i]]] + (D1[i]/
A[i])*(ArcTan[D1[i], (2*A[i] + B[i])] -
ArcTan[D1[i], B[i]]))];
If[ 4*A[i]*E1[i] - B[i]^2 == 0, PF2[i] = 0,
PF2[i] =
lg[i]*(nx[i]*(xb[i] - xi) + ny[i]*(yb[i] - eta))/
D1[i]*(ArcTan[D1[i], 2*A[i] + B[i]] - ArcTan[D1[i], B[i]])];
sum = sum + phi[i]*PF2[i] - dphi[i]*PF1[i];, {i, 1, NN}];
{xi, eta, sum/Pi}, {xi, .1, .9, .1}, {eta, 0.1, .9, .1}];


Equivalent FEM Code

reg = Rectangle[{0, 0}, {1, 1}];
bc = {DirichletCondition[u[x, y] == Cos[Pi y], x == 1],
DirichletCondition[u[x, y] == 0, x == 0]};
bc1 = NeumannValue[0, y == 0 || y == 1] ;
U = NDSolveValue[{-Laplacian[u[x, y], {x, y}] == bc1, bc},
u, {x, y} \[Element] reg]

sol1 = Table[{x, y, U[x, y]}, {x, 0.1, .9, .1}, {y, 0.1, .9, .1}];


Comparison of two solutions and the solution itself

dsol = Flatten[sol - sol1, 1][[All, 3]];
{ListPlot[dsol],
ContourPlot[U[x, y], {x, y} \[Element] reg, PlotLegends -> Automatic,
ColorFunction -> "Rainbow", PlotRange -> All, Contours -> 20]}


• Dear Alex great job well done, can we evaluate the relative errors to see FEM or BEM is better with less error? Mar 19, 2020 at 20:54
• There is no point in deciding which method is more accurate. This is just an example of BEM for pedagogical purposes, which I worked on for a couple of hours. Whereas Mathematica FEM with NDSolve is a high-tech product that took several years. On this simple problem, we see that the solutions differ by a maximum of -0.000872012 at sol=-0.505865 and sol1=-0.504993. Mar 20, 2020 at 12:31
• Thank you Alex, some researchers said BEM better than FEM, that why I was asked to evaluate the error.....Best regards Mar 20, 2020 at 12:46
• For 3D problems, BEM is more accurate than FEM with equal split (I tested this). But on 2d there can be equal opportunities, as in this example. The accuracy of the BEM also depends on the accuracy of the integration method. Mar 20, 2020 at 14:03
• Many thanks Alex. Mar 20, 2020 at 16:59