How to compute the inner integral with quad?

56 views Asked by At

I am reading the paper by Floris​ and try to compute the numerical value of the normalization constant C_I from eq. (20).

From the eqution

2]

and property of transition probability density function (PDF) enter image description here, the normalization constant is:

enter image description here

In my case, m_1(x)=-(a*x+b*x**3), m_2(x)= c+ sigma**2*x, a=1.0, b=0.5, c = 1.0, sigma = 1.0.

I have used the quad() function in my attempt for the inner integral:

import numpy as np
exp = np.exp
inf = np.inf
log  = np.log
sqrt = np.sqrt

from scipy.integrate import quad

a = 1.0;  b = 0.5; c = 1.0; sigma = 1.0

tspan = np.linspace(0.0, 2.5, 5001)
x0 = 0.1

def m_1(x, t):
  return -(a*x + b*x**3)

def m_2(x, t):
  return (c + sigma**2 * x)

def inner_integrand(x0, x, t):
  return 2 * m_1(x,t) / m_2(x,t)**2

inner = quad(inner_integrand, -inf, inf, args=(x0, tspan))[0]
print('inner_integrand=', inner)
# inner_integrand= -1.693896469019819

and I have obtained the result inner_integrand= 0.33223140495839776.

When I have used the wolfram alpha, the result is 'integral does not convergence'

Question. How to compute the inner integral?

1

There are 1 answers

1
julaine On

integrand will be called with some x selected by the algorithm and your extra arguments after that.

Your call quad(integrand, -inf, inf, args=(x0, tspan)) will make calls like integrand(x, x0, tspan) and while I do not know what integrand is, this does not seem to be what the function expects.