2

I faced an issue when trying to use quad to integrate a function. Essentially, I have two versions of code where I define t(a) in different places. Both codes looks the same to me but the result I am getting is slightly different.

I am guessing that it is due to the error associated with using the quad method, but am not too sure. Would appreciate any help!

import numpy as np
from scipy.integrate import quad

s = 0.05

# Version 1
def POi1(w):
    def t(a):
        return (1/(0.27*a**(-1)+(1-0.27)*a**(-3*w-1))**(1/2))
    return (np.exp(-((1+w)**2)/(2*s**2))*(1/(quad(t, 0, 1)[0]+((2/3)*(1/(np.abs(1+w)*1*(1-0.27)**2))))))

PO1 = quad(POi1, -np.inf, -1)[0]
print(PO1)

#Version 2
def t(a):
    return (1/(0.27*a**(-1)+(1-0.27)*a**(-3*w-1))**(1/2))

def POi1(w):
    return (np.exp(-((1+w)**2)/(2*s**2))*(1/(quad(t, 0, 1)[0]+((2/3)*(1/(np.abs(1+w)*1*(1-0.27)**2))))))

PO1 = quad(POi1, -np.inf, -1)[0]
print(PO1)
kederrac
  • 16,819
  • 6
  • 32
  • 55
Azuyugi
  • 23
  • 3
  • even if you set s=1 (any number) version 2 will not work since there is no 'w' variable for t(a) defined – kederrac Sep 12 '20 at 09:46
  • Hi kerderrac, thank you for your comments. I have edited my code. Made a mistake when copying over the code. If you try to run the edited code, you should see that the answers will be different. – Azuyugi Sep 12 '20 at 09:48
  • trying to run your code, same issue: NameError: name 'w' is not defined – kederrac Sep 12 '20 at 09:49
  • For version 2, the w will be defined in POi1 when I integrate t(a) using quad. Hopefully, this is clearer! Sorry if it wasn't earlier – Azuyugi Sep 12 '20 at 09:49
  • V1, PO1 = 0.0019037562389785248 V2, PO1 = 0.0019343431215966553 – Azuyugi Sep 12 '20 at 09:53
  • if your run on jupyter just restart your kernel and run only your code posted here, you will see that will not work, maybe you have some variable w already defiend – kederrac Sep 12 '20 at 09:55

1 Answers1

0

your current code doesn't work (version 2) because the function t will take as parameter a, only one parameter, but you can see that requires 2 params, a and w, it may work if you have a variable w already defined (maybe you work in jupyter notebook) and this it may be the cause of your different results

to make your code work and have the same result for your version 2 you can use:

def get_t(w):
    def t(a):
        return (1/(0.27*a**(-1)+(1-0.27)*a**(-3*w-1))**(1/2))
    
    return t

def POi1(w):
    t = get_t(w)
    return (np.exp(-((1+w)**2)/(2*s**2))*(1/(quad(t, 0, 1)[0]+((2/3)*(1/(np.abs(1+w)*1*(1-0.27)**2))))))

PO1 = quad(POi1, -np.inf, -1)[0]
print(PO1)

in this example a closure it is used to pass the w variable, you will get the same results

kederrac
  • 16,819
  • 6
  • 32
  • 55