CODEWITHSUNDEEP
pythonmathscipy
tux007
tux007

Reputation: 312

scipy.integrate.quad with non plausible values

He guys,

I'm al little bit in trouble with math. So there is an old question and it's not really solved. I thought about editing the old one, but I think it's even good to start a new question.

As you can see below there is an example of my experimental code. There are the functions kf and kfm. kf is a little bit funcy fancy math stuff from an engineer. kfm is just the integral of kf.

When I use the trapezoid rule, the value of kfm(0.05) is ~67. It's easy, because kf(0.05) is ~135 and the first step of the trapezoid is 1/2 of kf.

So I don't know why kfm (the integral of kf) is starting at zero? I think it is wrong!

'''python code '''

import matplotlib.pyplot as plt
from scipy import integrate
import math

# umf is just the x axis for kf and kfm
list_umf = [0.05, 0.1, 0.15000000000000002, 0.2, 0.25, 0.3, 0.35000000000000003, 0.4, 0.45, 0.5, 0.55,
            0.6000000000000001, 0.6500000000000001, 0.7000000000000001, 0.7500000000000001, 0.8]

# the real kf function
def kf(phi):
    k = 841.17
    m1 = -0.00112
    m2 = 0.17546
    m3 = 0.01271
    m4 = 0.00116
    phi_p = 3.42
    theta_einlauf = 1200
    if phi < 0.02:
        phi = 0.02
    return k * math.exp(m1 * theta_einlauf) * (phi ** m2) * (phi_p ** m3) * math.exp(m4 / phi)
list_kf_func = [kf(x) for x in list_umf]

# the real kfm function
def kfm(u, o):
    return integrate.quad(kf, u, o)[0]

u = list_umf[0]
list_kfm_func = [kfm(u, o) for o in list_umf]

# here are some plots
plt.plot(list_umf, list_kf_func, label='kf_func')
plt.plot(list_umf, list_kfm_func, label='kfm_func')
plt.legend(loc="upper left")

Here is the plot, generated with the code: Here is the plot, generated with the code

Upvotes: 1

Views: 121

Answers (1)

adrianop01
adrianop01

Reputation: 600

Your trapezoidal rule estimation is incorrect.

For the first scipy.integrate from lower limit a=0 to upper limit b=0.05, Using trapezoidal rule, enter image description here

f is but kf(). Thus calculating by trapezoidal rule: 0.5*(kf(0.05)+kf(0))*(0.05-0) gives 6.342386940782528.

In fact, if you check the output from your Python code, kfm(0,0.05) = list_kfm_func[0] = integrate.quad(kf, 0, 0.05)[0] = 6.204186103148345, which agrees with the trapezoidal rule results.

Edit

I saw your edit history of adding u = list_umf[0] in your question, i.e. lower limit a= 0.05. My old answer assumes a lower limit a = 0 as this code was missing in your original code, my apologies.

However, in this case, a = b = 0.05 for the first term of list_kfm_func, i.e. both numerical integration method shall give 0 as lower limit equals to upper limit, therefore, both integration still agrees with each other and your trapezoidal rule estimation stated in your question is still not correct.

By trapezoidal rule: 0.5*(kf(0.05)+kf(0.05))*(0.05-0.05) = 0

By scipy.integrate: kfm(0.05,0.05) = list_kfm_func[0] = integrate.quad(kf, 0.05, 0.05)[0] = 0.0

Upvotes: 3

Related Questions