# On the Dirichlet integral

Continuing this topic, once the Dini-Neumann problem is solved:

$$\begin{cases} \Delta f = 0 & \text{on} \; \Omega \\ \frac{\text{d}}{\text{d} n} f = y\,n_x - x\,n_y & \text{on} \; \partial\Omega \end{cases}$$

with $$n = (n_x,\,n_y)$$ exterior normal vector to $$\partial\Omega$$, the goal is the calculation of:

$$\color{blue}{J \equiv \iint_{\Omega} \left(f_x\,y - f_y\,x\right)\text{d}x\,\text{d}y}\,.$$

On the other hand, we also have:

$$J = \iint_{\Omega} \left[\left(f\,y\right)_x - \left(f\,x\right)_y\right]\text{d}x\,\text{d}y$$

which by the divergence theorem is equivalent to:

$$J = \int_{\partial\Omega} f\left(y\,n_x - x\,n_y\right)\text{d}s$$

and for the boundary condition:

$$J = \int_{\partial\Omega} f\left(\frac{\text{d}}{\text{d}n}f\right)\text{d}s = \int_{\partial\Omega} \left(f\,f_x\,n_x + f\,f_y\,n_y\right)\text{d}s\,.$$

This done, again by the divergence theorem:

$$J = \iint_{\Omega} \left[(f\,f_x)_x + (f\,f_y)_y\right]\text{d}x\,\text{d}y = \iint_{\Omega} \left(f_x^2 + f_y^2 + f\,\Delta f\right)\text{d}x\,\text{d}y$$

that is, since $$f$$ is an harmonic function on $$\Omega$$:

$$\color{\red}{J = \iint_{\Omega} \left(f_x^2 + f_y^2\right)\text{d}x\,\text{d}y}\,,$$

known in the literature as the Dirichlet integral.

Therefore, after determining the harmonic function:

Needs["NDSolveFEM"];

Ω = ToElementMesh[ImplicitRegion[x^4 + y^4 <= 1, {x, y}], "IncludePoints" -> {{0, 0}}];

f = NDSolveValue[{Laplacian[g[x, y], {x, y}] == NeumannValue[y (4 x^3) - x (4 y^3), True],
DirichletCondition[g[x, y] == 0, x == 0 && y == 0]}, g, {x, y} ∈ Ω];


and verified that it's accurate enough:

NIntegrate[f[x, y], {x, y} ∈ Ω]
NIntegrate[Laplacian[f[x, y], {x, y}], {x, y} ∈ Ω]


-0.00000293157

0.00292529

I tried to calculate the Dirichlet integral in its two equivalent versions:

NIntegrate[D[f[x, y], x] y - D[f[x, y], y] x, {x, y} ∈ Ω]
NIntegrate[D[f[x, y], x]^2 + D[f[x, y], y]^2, {x, y} ∈ Ω]


-0.48313

1.84564

but it's evident that something went wrong, but I don't know what! Ideas? Thank you!

EDIT: Thanks to the precious comment of xzczd I understood that:

• $$(n_x,n_y)$$ must necessarily have unitary length;

• before the Laplace operator it must absolutely be a minus sign.

In fact, writing:

c[x_, y_] = x^4 + y^4;
{nx, ny} = Normalize[Grad[c[x, y], {x, y}]];

Needs["NDSolveFEM"];
Ω = ToElementMesh[ImplicitRegion[c[x, y] <= 1, {x, y}],
"IncludePoints" -> {{0, 0}}];

f = NDSolveValue[{-Laplacian[g[x, y], {x, y}] == NeumannValue[y nx - x ny,
True], DirichletCondition[g[x, y] == 0, x == 0 && y == 0]},
g, {x, y} ∈ Ω];

NIntegrate[f[x, y], {x, y} ∈ Ω]
NIntegrate[Laplacian[f[x, y], {x, y}], {x, y} ∈ Ω]

NIntegrate[D[f[x, y], x] y - D[f[x, y], y] x, {x, y} ∈ Ω]
NIntegrate[D[f[x, y], x]^2 + D[f[x, y], y]^2, {x, y} ∈ Ω]


we get:

0.000000504282

-0.000799624

0.126589

0.126595

which is approximately what I want. The only thing that would remain to understand is whether it were possible to obtain a better approximation of those integrals. Thanks again!

• Can you add a reference for Dirichlet integral? It doesn't seem to belong to those mentioned in wiki page: en.wikipedia.org/wiki/Dirichlet_integral Dec 5, 2022 at 2:07
• For divergence theorem, $(n_x,n_y)$ should be unit vector, isn't it? Dec 5, 2022 at 8:34

Let's make the mesh denser by setting MaxCellMeasure option:

Ω =
ToElementMesh[ImplicitRegion[c[x, y] <= 1, {x, y}], "IncludePoints" -> {{0, 0}},
MaxCellMeasure -> 0.001];


Then the output of

NIntegrate[D[f[x, y], x] y - D[f[x, y], y] x, {x, y} ∈ Ω]
NIntegrate[D[f[x, y], x]^2 + D[f[x, y], y]^2, {x, y} ∈ Ω]


is

0.126603
0.126603