I have a function:

delt=1 #trial
def f(z):
    return ((1-2*z)*np.exp(-delt/z))/(((1-z)**(2+delt))*(z**(2-delt)))

I also have a variable:

import scipy.integrate as integrate
var=integrate.quad(f,0,0.5)[0]  # equals 0.040353419593637516

Now I am trying to find the value p such that

integrate.quad(f,0.5,p)= var

and manually I can check that it is around 0.605

I define the following function to be used in optimization:

def integral(p):
    return integrate.quad(f,0.5, p)[0]-var

However, i get the following results:

import scipy.optimize as op
In[26]: op.root(integral,0.61)
    fjac: array([[-1.]])
     fun: -0.040353420516861596
 message: 'The iteration is not making good progress, as measured by the \n  improvement from the last ten iterations.'
    nfev: 18
     qtf: array([ 0.04035342])
       r: array([ 0.00072888])
  status: 5
 success: False
       x: array([ 0.50002065])
In[27]: op.fsolve(integral,0.61)
Out[27]: array([ 0.50002065])

Any idea why both root and fsolve might be failing?

  • 1
    $\begingroup$ integrate.quad(f, 0.5, 0.605) is approximately 0 and not close to var as you claim. Also, a plot of integral(p) for p>0.5 suggests that it is less than zero $\endgroup$
    – Stelios
    Mar 9 '17 at 9:06

As @Stelios mentioned, It seems like the integral from 0.5 to ~0.605 is close to 0, and then it turns negative.The antiderivative of your function is given by

$$\int f(z; t=1)\, dz = \frac{e^{-\frac{1}{z} - 1}\left[e\, (1-2z)z - e^{1/z}(z-1)^2\operatorname{Ei}\left(\frac{z-1}{z}\right)\right]}{2(z-1)^2}\, ,$$

keep in mind that it does not work for $z=0$, indeed, your original function is not defined at $z=0$, and the limit does not even exist. If you make the plot of the antiderivative you can see these features

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.