Wrong integral approximation for a singular function

dolfinx
1

I have a problem when I simply try to integrate the function 1/x in a simple way. I have tried to change the degree of squareness of the method but the result is still very bad. The same happens when I refine the mesh.
Is there another way to integrate correctly or am I doing something wrong?
I attach a small code with an example and an image showing the variation of the result in relation to the degree of squareness and the number of elements.
Thanks in advance.

    Nx = 6  # number of elements
    # Meshing the defined domain using a Nx nodes resolution
    mesh = msh.create_interval(MPI.COMM_WORLD, Nx, [0, 1.0])
    dx = ufl.Measure("dx", domain=mesh, metadata={"quadrature_degree": 10},)
    x = ufl.SpatialCoordinate(mesh)
    V = femx.FunctionSpace(mesh, ("CG", 1))
    v = ufl.TestFunction(V)
    gamma_s = lambda x: 1/x[0]
    l = femx.form(gamma_s(x) * ufl.conj(v) * dx)
    b = femx.petsc.assemble_vector(l)
    A = np.sum(b.array[:])
2

Isn’t that integral \int_0^1 \frac{dx}{x}? That’s undefined when it has the lower limit at 0, it doesn’t converge.

3

Hi @dparsons, thank you for your answer. Yes, this is the integral. As you said it’s undefined at 0, but its integral should increase exponentially as long as the accuracy increases, isn’t? Exploding beahviour should be observed.
That’s the picture I want to share in my previous message:
test_quadrature_degree_singular_function

4

I see what you mean. It is increasing though. I can’t interpret the details of the numerics myself, perhaps others can comment.

5

I think logarithmic (rather than exponential) growth of the integral is to be expected. Assume that the quadrature scheme can exactly compute the integral everywhere except in the first element (near the singularity at x=0). Then we can write the numerical approximation of the integral, I, as the following (where L_e is the length of an element):

I = \int_0^1{\frac{dx}{dx}} = \int_0^{L_e} {\frac{dx}{x}} + \int_{L_e}^1 {\frac{dx}{x}}

The integral over the first element can be approximated numerically by transforming to an integral over the interval [-1, 1] (with x = L_e (\xi + 1)/2) and replacing with a sum over the quadrature points \xi_i and weights w_i:

\int_0^{L_e} {\frac{dx}{x}} = \int_{-1}^1 {\frac{d\xi}{\xi + 1}} \approx \sum_{i=1}^n {\frac{w_i}{\xi_i + 1}} = A

Note that this is a constant with no dependence on the mesh size. The remaining term of the integral can be integrated exactly:

\int_{L_e}^1 {\frac{dx}{x}} = \ln{1} - \ln{L_e} = - \ln{L_e}

Finally, making the substitution L_e = 1/N_e, where N_e is the number of elements, we obtain:

I = A + \ln{N_e}\,,

i.e. the approximation of the integral should grow logarithmically with the number of elements. Loosely speaking, this is much slower than a linear growth, and much much (= much²) slower than an exponential growth. This appears to match your results; running

import numpy as np
import ufl
from dolfinx import fem, mesh
from mpi4py import MPI
from petsc4py import PETSc

q = 10
Nx = 64
for Nx in range(10):
    Ne = 2**Nx
    msh = mesh.create_interval(MPI.COMM_WORLD, Ne, [0, 1.0])
    dx = ufl.Measure("dx", domain=msh, metadata={"quadrature_degree": q},)
    x = ufl.SpatialCoordinate(msh)
    V = fem.FunctionSpace(msh, ("CG", 1))
    v = ufl.TestFunction(V)
    l = fem.form(1/x[0] * ufl.conj(v) * dx)
    b = np.sum(fem.petsc.assemble_vector(l).array[:])
    print(Ne, b)

produces

1 4.899999999999997
2 5.593147179886525
4 6.286294360445773
8 6.979441541005717
16 7.672588721565662
32 8.365735902125607
64 9.05888308268555
128 9.752030263245496
256 10.445177443805441
512 11.138324624365389

and plotting the results

import numpy as np
from matplotlib import pyplot as plt
data=np.array([[1,4.899999999999997],
[2,5.593147179886525],
[4,6.286294360445773],
[8,6.979441541005717],
[16,7.672588721565662],
[32,8.365735902125607],
[64,9.05888308268555],
[128,9.752030263245496],
[256,10.445177443805441],
[512,11.138324624365389]])
plt.plot(np.log(data[:,0]), data[:,1])
plt.xlabel("$\\log(N_{elem})$")
plt.ylabel("$\\int_0^1{dx/x}$")
plt.show()

gives
Figure_1

(The accuracy of the quadrature scheme outside the first element can be verified with the following code:)

import numpy as np
import ufl
from dolfinx import fem, mesh
from mpi4py import MPI
from petsc4py import PETSc

q = 10
Nx = 64
for Nx in range(1, 10):
    Ne = 2**Nx
    msh = mesh.create_interval(MPI.COMM_WORLD, Ne-1, [1/Ne, 1.0])
    dx = ufl.Measure("dx", domain=msh, metadata={"quadrature_degree": q},)
    x = ufl.SpatialCoordinate(msh)
    V = fem.FunctionSpace(msh, ("CG", 1))
    v = ufl.TestFunction(V)
    l = fem.form(1/x[0] * ufl.conj(v) * dx)
    b = np.sum(fem.petsc.assemble_vector(l).array[:])
    print(Ne, b, np.log(Ne))
3 Likes
6

Thank you @conpierce8 . This definitely explains the results I had obtained. I think that as I was told in Slack custom definition for ufl.dx may let me obtain better results.
I share the answers I had obtained: