I am trying to solve a simple integral in Python using different methods for numerical integration. However, when I try to calculate the integrals, I obtain wrong results for small values of my x coordinate. This is a messy because I am trying to solve some system of equations in which one of the variables I am evolving on time is the integral of some quantities. Of course, I don’t have a functional form for the expression I want to integrate but instead I have an array of data (one f[i] for each x[i]). Do you know a way to obtain more accurate results of my integrals?

Here, I present a simple example of the problem I am having. Suppose I have some data that, for now, I will only consider a simple x^7 (in general I don’t have a functional form for the functions I want to integrate, with this is just an example) and suppose I want to integrate such a function. The analytic solution of the integral is simple Int = x^8/8. Then, if I print or plot 8*Int/x^8 what I should have is that each entry of the output should be equal to one. Then, let me consider the following code:

import numpy as np
from scipy.integrate import simps 

x = np.linspace(0,10,101)
f = x**7
Int = np.array([simps(f[:i+1], x[:i+1]) for i in range(len(x))])
Div = 8*Int/(x**8)
Div[0] = 1 #To avoid division over zero
print(Div)

The output of this simple code is:

[1.           4.           1.375        1.2687090383 1.033203125
1.05997312   1.0069158665 1.0195120352 1.0022277832 1.0080582909
1.00092      1.0038955372 1.0004456376 1.0021034429 1.0002411835
 1.0012324471 1.0001416206 1.0007688345 1.0000885172 1.0005040075
 1.000058125  1.0003439529 1.0000397249 1.0002426572 1.0000280617
 1.0001760417 1.000020381  1.0001307892 1.000015157  1.0000991831
 1.0000115043 1.000076571  1.0000088885 1.0000600497 1.0000069756
 1.0000477527 1.0000055507 1.0000384478 1.0000044717 1.0000313025
 1.0000036426 1.0000257424 1.000002997  1.0000213637 1.0000024883
 1.0000178775 1.0000020831 1.0000150742 1.0000017571 1.0000127994
 1.0000014925 1.0000109379 1.0000012758 1.0000094029 1.0000010971
 1.0000081279 1.0000009486 1.0000070618 1.0000008244 1.000006165
 1.0000007199 1.000005406  1.0000006314 1.0000047604 1.0000005561
 1.0000042083 1.0000004917 1.000003734  1.0000004364 1.0000033246
 1.0000003886 1.0000029699 1.0000003472 1.0000026612 1.0000003112
 1.0000023916 1.0000002797 1.0000021553 1.0000002521 1.0000019475
 1.0000002278 1.0000017641 1.0000002064 1.0000016018 1.0000001874
 1.0000014577 1.0000001706 1.0000013295 1.0000001556 1.0000012151
 1.0000001422 1.0000011127 1.0000001303 1.0000010209 1.0000001195
 1.0000009384 1.0000001099 1.000000864  1.0000001012 1.0000007969
 1.0000000933]
 

As you can see, most of the results are very close to 1 (as expected) except for the points that are very close to the point x = 0. I need to have a good control of these integrals because I am using them for solving numerically a system of equations. The problem of having these bad results is that in my numerical scheme it ends giving me some gradients that ends destroying my numerical solutions. So, do you know a way to deal with this?

PS: If instead of using simps I use mpmath, I can obtain the correct result. However, mpmath works only for functions that we already know (as it is the case in this simple example of a x^7 function). However, as I previously mentioned, in reality I don’t know a functional form of the expression I wan’t to integrate and instead I only have an array of data.

PS2: I can also try to increase the size of my array x. But even if I do this, I always have errors in the first points closest to x = 0, and those errors end up ruining my numerical solutions.

Best regards, LP