# Solving the Liouville-Bratu-Gelfand equation on the surface of a cube

How do I solve the equation:

Laplacian[ u[x,y,z], {x,y,z} ] = (Exp[ u[x,y,z] + h[x,y,z] ]) -1, where h[x,y,z]= Sum [ { Log |{x,y,z} - {x_i, y_i, z_i}|^2 }, { i, 1, N}], on the surface of a cube, which is a compact surface without a boundary?

• Did the edit I made to this post not solve this already? – Henrik Schumacher Dec 27 '18 at 19:41
• I didn't notice the edit you made before. I thought you were still working on it. I changed the function h in the edit you made. I'm running the code to see if it works with these edits. Please see the edit I made in the previous post to see the edited code I'm running. – Thando Dec 27 '18 at 20:22
• Is the Sum supposed to range over only a few points or over all points on the surface? Actually, I am missing quite a lot detail here. – Henrik Schumacher Dec 27 '18 at 20:34
• For a numerical model, add the number of points and the points xi,yi themselves, eg random points. – Alex Trounev Dec 27 '18 at 20:45
• @Henrik, the Sum ranges over random points on the surface. Each of these random points represents a vortex position and there are N of these. So N ranges from 1 (which is one vortex) to any integer which of course must be <= the total number of points on the surface of the cube. – Thando Dec 27 '18 at 20:52

I will give an explanation to the solution of the problem. First, $$u[x,y]=-h[x,y]$$ is the exact solution of the problem on a plane for h[x_, y_] := Sum[Log[(x - xn[[i]])^2 + (y - yn[[i]])^2], {i, 1, n}]. This raises the question of how the influence of the source given on one face extends to all other faces? If it is through 3D, then this is one problem, and if through 2D, then this is another problem and it has a different solution. Secondly, it is possible to cut and unfold the surface of the cube on a plane. Therefore, the problem is reduced to 2D. But it is necessary to sew the solution on all cuts. Third, we can solve the problem on the surface of a cube in 3D, as did Henrik Schumacher in this post. I wrote code for 2D sources that are a solution to the Laplace equation in 2D. I did not continue the sources through the cuts. If the author needs a continuation, I will add a couple of lines. It shows what the solution looks like on a plane.

c1 = ImplicitRegion[-1 <= x <= 0 && 0 <= y <= 1, {x, y}];
c2 = ImplicitRegion[0 <= x <= 1 && 0 <= y <= 1, {x, y}];
c3 = ImplicitRegion[1 <= x <= 2 && 0 <= y <= 1, {x, y}];
c4 = ImplicitRegion[2 <= x <= 3 && 0 <= y <= 1, {x, y}];
c5 = ImplicitRegion[0 <= x <= 1 && 1 <= y <= 2, {x, y}];
c6 = ImplicitRegion[0 <= x <= 1 && -1 <= y <= 0, {x, y}];
cubS = RegionUnion[c1, c2, c3, c4, c5, c6];

n = 20; n1 = RandomInteger[{0, n}]; n2 =
RandomInteger[{0, n - n1}]; n3 =
RandomInteger[{0, n - n2 - n1}]; n4 =
RandomInteger[{0, n - n3 - n2 - n1}]; n5 =
RandomInteger[{0, n - n4 - n3 - n2 - n1}];
n6 = n - n5 - n4 - n3 - n2 - n1;
t0 = 1/50; k = 10; x1 = Table[RandomReal[{-1, 0}], {n1}]; y1 =
Table[RandomReal[{0, 1}], {n1}]; x2 =
Table[RandomReal[{0, 1}], {n2}]; y2 = Table[RandomReal[{0, 1}], {n2}];
x3 = Table[RandomReal[{1, 2}], {n3}]; y3 =
Table[RandomReal[{0, 1}], {n3}];
x4 = Table[RandomReal[{2, 3}], {n4}]; y4 =
Table[RandomReal[{0, 1}], {n4}];
x5 = Table[RandomReal[{0, 1}], {n5}]; y5 =
Table[RandomReal[{1, 2}], {n5}];
x6 = Table[RandomReal[{0, 1}], {n6}]; y6 =
Table[RandomReal[{-1, 0}], {n6}];
xn = Join[x1, x2, x3, x4, x5, x6]; yn =
Join[y1, y2, y3, y4, y5, y6]; pn =
Table[{xn[[i]], yn[[i]]}, {i, 1, Length[xn]}];
Show[RegionPlot[cubS, AspectRatio -> Automatic],
Graphics[{Red, PointSize[.01], Point[pn]}]]

h[x_, y_] := Sum[Log[(x - xn[[i]])^2 + (y - yn[[i]])^2], {i, 1, n}]

U[0][x_, y_] := -h[x, y]
Do[U[i] =
NDSolveValue[{(u[x, y] - U[i - 1][x, y])/t0 ==
Laplacian[u[x, y], {x, y}] - Exp[U[i - 1][x, y] + h[x, y]] + 1,
DirichletCondition[u[x, y] == U[i - 1][-1, y],
x == 3 && 0 <= y <= 1],
DirichletCondition[u[x, y] == U[i - 1][0, x],
y == 0 && -1 <= x <= 0],
DirichletCondition[u[x, y] == U[i - 1][1, -x + 1],
y == 0 && 1 <= x <= 2],
DirichletCondition[u[x, y] == U[i - 1][0, -x + 1],
y == 1 && -1 <= x <= 0],
DirichletCondition[u[x, y] == U[i - 1][1, x],
y == 1 && 1 <= x <= 2],
DirichletCondition[u[x, y] == U[i - 1][x - 2, -1],
y == 0 && 2 <= x <= 3],
DirichletCondition[u[x, y] == U[i - 1][x - 2, 2],
y == 1 && 2 <= x <= 3]}, u, {x, y} \[Element] cubS];, {i, 1,
k}];
Table[ContourPlot[U[i][x, y], {x, y} \[Element] cubS, Contours -> 20,
ColorFunction -> Hue, PlotLegends -> Automatic, PlotLabel -> i,
PlotRange -> All, AspectRatio -> Automatic], {i, 1, k}]


The distribution of sources on the surface and solution evolution starting from i= 0 and up toi = 10. Note that the solution $$u(x,y)=-h(x,y)$$ at i=0 contains logarithmic singularities.

• I see what you mean Alex. Your solution is equivalent to Henrik's. This is great! It makes things clearer. Thank you for that. – Thando Dec 30 '18 at 8:23
• This is similar to a 3D solution due to splicing through cuts, but there is a difference. First, the sources in 2D are singular (these are vortices). Therefore, the numerical solution of the equation in 2D is unstable. Secondly, the influence of sources in 2D is different from 3D. I will do the calculations in 2D and in 3D and show the differences. – Alex Trounev Dec 30 '18 at 11:14
• Thanks Alex. I'm looking forward to seeing your results. The comparison is very interesting. – Thando Dec 30 '18 at 15:54