I receive three different warnings in scipy integration in python

Question

I have the code structure below. I would like to get some numerical results here.

import numpy as np
import scipy
from scipy import integrate

alpha = .99
t = np.linspace(0, .85, 5)
s = np.empty_like(t)
f = np.power(t - s, -alpha)
Int = integrate.simpson(f, s)
Int

I got the warnings below. I understand that the first term in t, that is t[0] causes the errors, particularly first two warnings. But I do not know how I can avoid these warnings. I cannot change the alpha,t or f.

<ipython-input-1-6b0d0757bfac>:8: RuntimeWarning: invalid value encountered in power
  f = np.power(t-s, -alpha)
/usr/local/lib/python3.8/dist-packages/scipy/integrate/_quadrature.py:414: RuntimeWarning: invalid value encountered in true_divide
  h0divh1 = h0 / h1
/usr/local/lib/python3.8/dist-packages/scipy/integrate/_quadrature.py:416: RuntimeWarning: invalid value encountered in true_divide
  y[slice1] * (hsum * hsum / hprod) +
nan

I tried to take t = np.linspace(1e-8, .85, 5). It did not work.

EDIT: The t is stable. I cannot change it. I need to derive s from t or find a new definition for s. So, I have

t = np.linspace(0, .85, 5)
array([0.    , 0.2125, 0.425 , 0.6375, 0.85  ])

Let us take the s as an array of zeros (array of ones also does not work for s)

s = np.zeros_like(t)
array([0., 0., 0., 0., 0.])

and add 1 to f to avoid the warning of RuntimeWarning: divide by zero encountered in power

f = (t1+1-s)** -.99
array([1.        , 0.82633295, 0.70424421, 0.61370619, 0.54387612])

After this when I used the simpson, it gives nan.

Int = integrate.simpson(f, s)
nan

Depending on the definition of s, I am constantly encountering a warning or an error .

The question is: What could be the right definition of s along with the t defined above?

Answer 1

There is an issue when you are using s because you created it like this:

s = np.empty_like(t)

What you are doing here is creating an array that has the same properties as t but is uninitialised, meaning that its shape and type are the same as t but its values are not set (thus "random").

For example, when I print s I have:

array([inf, nan, nan, nan, nan])

And you could have basically any values here.

So, you need to provide values to s before using it.

If you want to initialise s to zeros, you should replace that line by:

s = np.zeros_like(t)

Answer 2

Your first lines, and the results. With a size of 5, there's not excuse to not look at the arrays in full. No summary needed.

In [80]: alpha = .99
    ...: t = np.linspace(0, .85, 5)
    ...: s = np.empty_like(t)
    ...: f = np.power(t - s, -alpha)
C:\Users\paul\AppData\Local\Temp\ipykernel_2648\1323274646.py:4: RuntimeWarning: invalid value encountered in power
  f = np.power(t - s, -alpha)

5 equally space numbers:

In [81]: t
Out[81]: array([0.    , 0.2125, 0.425 , 0.6375, 0.85  ])

5 values as well; with empty they aren't predictable. sometimes I see 0s, but the is a good example of ridiculous "random" values. We don't want to use the s as is; it has to overwritten with something:

In [82]: s
Out[82]: 
array([6.36598737e-314, 1.90979621e-313, 8.48798317e-314, 1.69759663e-313,
       1.27319747e-313])

In [83]: f
Out[83]: array([       nan, 4.63355855, 2.33289375, 1.56158135, 1.17456015]

And the term that produced the warning (not error):

In [89]: (t[0]-s[0])
Out[89]: -6.365987374e-314
In [90]: (t[0]-s[0])**-.99
C:\Users\paul\AppData\Local\Temp\ipykernel_2648\147715877.py:1: RuntimeWarning: invalid value encountered in double_scalars
  (t[0]-s[0])**-.99
Out[90]: nan

Ekaba's suggestion:

In [91]: alpha = .99
    ...: t = np.linspace(0, .85, 5)
    ...: s = np.empty_like(t)
    ...: f = np.empty_like(t)
    ...: 
    ...: for i in range(1, len(t)):
    ...:     s[i] = t[i-1]
    ...:     f[i] = (t[i] - s[i])**(-alpha)
    ...:     

This overwrites s and f, but I'm not sure those are useful.

same 5 term t:

In [92]: t
Out[92]: array([0.    , 0.2125, 0.425 , 0.6375, 0.85  ])

s is now just t values, offset by 1. The first is still "random":

In [93]: s
Out[93]: 
array([-6.36598737e-314,  0.00000000e+000,  2.12500000e-001,
        4.25000000e-001,  6.37500000e-001])

And for f, again ignoring the first, the rest as the same:

In [94]: f
Out[94]: 
array([2.12199579e-314, 4.63355855e+000, 4.63355855e+000, 4.63355855e+000,
       4.63355855e+000])

The successive differences of t raised to -alpha:

In [96]: np.diff(t)
Out[96]: array([0.2125, 0.2125, 0.2125, 0.2125])    
In [97]: np.diff(t)**-.99
Out[97]: array([4.63355855, 4.63355855, 4.63355855, 4.63355855])

The integral is the area of a rectangle f[-1] high, and s[-1] wide:

In [104]: f[-1]*s[-1]
Out[104]: 2.953893574179393

I'm guessing the OP thought that he defined f as a function, that, simpson was going to pass successive values of s to it. Other integrate functions do take a function (e.g. quad), but this takes an array.

Simpson example

Here's an example that may be similar to what you want:

Specify an array of sample points:

In [109]: x = np.linspace(0,1,11)

and a t and f. Both are calculated from x. t adds a 1 to avoid the x=0 problem:

In [110]: t = (x+1)**2    
In [111]: f = (t-x)**-.99

Displaying the values:

In [112]: x,t,f
Out[112]: 
(array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ]),
 array([1.  , 1.21, 1.44, 1.69, 1.96, 2.25, 2.56, 2.89, 3.24, 3.61, 4.  ]),
 array([1.        , 0.90184157, 0.80818825, 0.72179746, 0.64388254,
        0.57463534, 0.51364905, 0.46021453, 0.41350815, 0.37270087,
        0.33701556]))

And the simpson integration:

In [113]: simpson(f,x)
Out[113]: 0.6073410199915571

Without plotting (or using sympy) I can only estimate whether that makes sense. The area under a curve that starts at 1 and decresses to .3, over the x range of 0 to 1 - .6 is a reasonable value.