Python: how to integrate functions with two unknown parameters numerically

Question

Now I have two functions respectively are

rho(u) = np.exp( (-2.0 / 0.2) * (u**0.2-1.0) )

psi( w(x-u) ) = (1/(4.0 * math.sqrt(np.pi))) * np.exp(- ((w * (x-u))**2) / 4.0) * (2.0 - (w * (x-u))**2)

And then I want to integrate 'rho(u) * psi( w(x-u) )' with respect to 'u'. So that the integral result can be one function with respect to 'w' and 'x'.

Here's my Python code snippet as I try to solve this integral.

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

x = np.linspace(0,10,1000)
w = np.linspace(0,10,500)

u = np.linspace(0,10,1000)

rho = np.exp((-2.0/0.2)*(u**0.2-1.0))

value = np.zeros((500,1000),dtype="float32")

# Integrate the products of rho with 
# (1/(4.0*math.sqrt(np.pi)))*np.exp(- ((w[i]*(x[j]-u))**2) / 4.0)*(2.0 - (w[i]*(x[j]-u))**2)
for i in range(len(w)):
    for j in range(len(x)):
        value[i,j] =value[i,j]+ integrate.simps(rho*(1/(4.0*math.sqrt(np.pi)))*np.exp(- ((w[i]*(x[j]-u))**2) / 4.0)*(2.0 - (w[i]*(x[j]-u))**2),u)

plt.imshow(value,origin='lower')
plt.colorbar()

As illustrated above, when I do the integration, I used nesting for loops. We all know that such a way is inefficient.

So I want to ask whether there are methods not using for loop.

Patol75 · Accepted Answer · 2020-10-28T01:26:21.980

2

Here is a possibility using scipy.integrate.quad_vec. It executes in 6 seconds on my machine, which I believe is acceptable. It is true, however, that I have used a step of 0.1 only for both x and w, but such a resolution seems to be a good compromise on a single core.

from functools import partial
import matplotlib.pyplot as plt
from numpy import empty, exp, linspace, pi, sqrt
from scipy.integrate import quad_vec
from time import perf_counter


def func(u, x):
    rho = exp(-10 * (u ** 0.2 - 1))
    var = w * (x - u)
    psi = exp(-var ** 2 / 4) * (2 - var ** 2) / 4 / sqrt(pi)
    return rho * psi


begin = perf_counter()
x = linspace(0, 10, 101)
w = linspace(0, 10, 101)
res = empty((x.size, w.size))
for i, xVal in enumerate(x):
    res[i], err = quad_vec(partial(func, x=xVal), 0, 10)
print(f'{perf_counter() - begin} s')
plt.contourf(w, x, res)
plt.colorbar()
plt.xlabel('w')
plt.ylabel('x')
plt.show()

UPDATE

I had not realised, but one can also work with a multi-dimensional array in quad_vec. The updated approach below enables to increase the resolution by a factor of 2 for both x and w, and keep an execution time of around 7 seconds. Additionally, no more visible for loops.

import matplotlib.pyplot as plt
from numpy import exp, mgrid, pi, sqrt
from scipy.integrate import quad_vec
from time import perf_counter


def func(u):
    rho = exp(-10 * (u ** 0.2 - 1))
    var = w * (x - u)
    psi = exp(-var ** 2 / 4) * (2 - var ** 2) / 4 / sqrt(pi)
    return rho * psi


begin = perf_counter()
x, w = mgrid[0:10:201j, 0:10:201j]
res, err = quad_vec(func, 0, 10)
print(f'{perf_counter() - begin} s')
plt.contourf(w, x, res)
plt.colorbar()
plt.xlabel('w')
plt.ylabel('x')
plt.show()

Addressing the comment

Just add the following lines before plt.show() to have both axes scale logarithmically.

plt.gca().set_xlim(0.05, 10)
plt.gca().set_ylim(0.05, 10)
plt.gca().set_xscale('log')
plt.gca().set_yscale('log')

edited Oct 28 '20 at 01:26

answered Oct 27 '20 at 07:30

Patol75

4,342
1
17
28

Thank you for your answer. By the way, is there a way to plot ```res``` w.r.t ```Log w``` and ```Log x```? – Stephen Wong Oct 28 '20 at 01:09
I have added some information with regard to having figure axes scale logarithmically. If, additionally, you want the axis to represent the logarithm of the variables rather than the variables, well, you can simply provide the logarithms in the call to `contourf`. – Patol75 Oct 28 '20 at 01:29
What does ```partial``` mean here? – Stephen Wong Oct 28 '20 at 02:56
You do not need it, but I was using it to change the value of `x` given to `func` at each iteration of the loop. [Documentation](https://docs.python.org/3/library/functools.html#functools.partial) | [StackOverflow question](https://stackoverflow.com/questions/15331726/how-does-functools-partial-do-what-it-does) – Patol75 Oct 28 '20 at 04:37
What is the base of ```log```, 'e' or '10'? – Stephen Wong Oct 28 '20 at 11:07
1

Please just have a look at the documentation for simple questions like this one. Extended discussions are also not supposed to be included in the comment section of an answer. Please open a new question if you need further help. – Patol75 Oct 28 '20 at 12:52
Why don't you input the argument ```u```, when you call the function ```func```? – Stephen Wong Oct 28 '20 at 13:40
1

It is happening under the hood when `quad_vec` executes. You would have noticed I still had to give the bounds of integration for `u`. Please refer to the documentation for further detail. – Patol75 Oct 28 '20 at 22:44

Python: how to integrate functions with two unknown parameters numerically

1 Answers1