# Wrong integral approximation for a singular function **URL:** \ **Category:** dolfinx\ **Created:** [September 19, 2022, 10:52am UTC](https://fenicsproject.discourse.group/t/wrong-integral-approximation-for-a-singular-function/9307 "2022-09-19T10:52:44Z")\ **Posts on this page:** 6\ **Page:** 1
**Author:** ![recon\_19](https://avatars.discourse-cdn.com/v4/letter/r/91b2a8/32.png) [@recon\_19](https://fenicsproject.discourse.group/u/recon_19)\ **Post date:** [September 19, 2022, 10:52am UTC](https://fenicsproject.discourse.group/t/wrong-integral-approximation-for-a-singular-function/9307/1 "2022-09-19T10:52:44Z")
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. ```auto 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[:]) ``` ---
**Author:** ![dparsons](https://avatars.discourse-cdn.com/v4/letter/d/5e9695/32.png) [@dparsons](https://fenicsproject.discourse.group/u/dparsons)\ **Post date:** [September 19, 2022, 11:58am UTC](https://fenicsproject.discourse.group/t/wrong-integral-approximation-for-a-singular-function/9307/2 "2022-09-19T11:58:07Z")
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. ---
**Author:** ![recon\_19](https://avatars.discourse-cdn.com/v4/letter/r/91b2a8/32.png) [@recon\_19](https://fenicsproject.discourse.group/u/recon_19)\ **Post date:** [September 19, 2022, 12:08pm UTC](https://fenicsproject.discourse.group/t/wrong-integral-approximation-for-a-singular-function/9307/3 "2022-09-19T12:08:39Z")
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](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/5/5ab1290999a41e9b99c92ef0df390de892706071.png) ---
**Author:** ![dparsons](https://avatars.discourse-cdn.com/v4/letter/d/5e9695/32.png) [@dparsons](https://fenicsproject.discourse.group/u/dparsons)\ **Post date:** [September 20, 2022, 9:09am UTC](https://fenicsproject.discourse.group/t/wrong-integral-approximation-for-a-singular-function/9307/4 "2022-09-20T09:09:48Z")
I see what you mean. It is increasing though. I can’t interpret the details of the numerics myself, perhaps others can comment. ---
**Author:** ![conpierce8](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/conpierce8/32/2464_2.png) [@conpierce8](https://fenicsproject.discourse.group/u/conpierce8)\ **Post date:** [September 20, 2022, 5:08pm UTC](https://fenicsproject.discourse.group/t/wrong-integral-approximation-for-a-singular-function/9307/5 "2022-09-20T17:08:12Z")
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 ```python 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 ```nohighlight 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 ```python 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](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/7/7b29fd9a2126d070d7f931bc9ced3363510d9736.png) (The accuracy of the quadrature scheme outside the first element can be verified with the following code:) ```python 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)) ``` ---
**Author:** ![recon\_19](https://avatars.discourse-cdn.com/v4/letter/r/91b2a8/32.png) [@recon\_19](https://fenicsproject.discourse.group/u/recon_19)\ **Post date:** [September 22, 2022, 7:09am UTC](https://fenicsproject.discourse.group/t/wrong-integral-approximation-for-a-singular-function/9307/6 "2022-09-22T07:09:41Z")
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: > **[Slack is your digital HQ](https://fenicsproject.slack.com/?redir=%2Farchives%2FC9EQQU29M%2Fp1663594486921949)** > > Slack is a new way to communicate with your team. It’s faster, better organized, and more secure than email.