In order to solve a equation where the left hand side appears under an integral on the right hand side:
$$ B(p^2) = C\int_0^{p^2}f_1\left(B(q^2),q^2\right)\mathrm{d}q^2 + C\int_{p^2}^{\Lambda^2} f_2\left(B(q^2),q^2\right)\mathrm{d}q^2 $$
I have written the following code to do solve this equation iteratively:
import numpy as np
import time
def solver(alpha = 2.):
L2 = 1 # UV Cut-off
N = 1000 # Number of grid points
maxIter = 100 # Maximum no. of iterations
# The grid on which the function is evaluated
p2grid = np.logspace(-10, 0, N, base=10)
dq = np.array([0] + [p2grid[i+1]-p2grid[i] for i in range(len(p2grid)-1)])
Bgrid = np.ones(N)
Bgrid_new = np.empty(N) # Buffer variable for new values of B
C = alpha / (4*np.pi)
for j in range(maxIter):
for i, p2 in enumerate(p2grid):
# If lower and upper limit of an integral are the same, set to 0
if i > 0:
n = i + 1
int1 = np.add.reduce((p2grid[0:n] * 3*Bgrid[0:n] /\
(p2 * (p2grid[0:n] + Bgrid[0:n]**2))) * dq[0:n])
else:
int1 = 0
if i < N - 1:
int2 = np.add.reduce((3*Bgrid[i:] /
(p2grid[i:] + Bgrid[i:]**2)) * dq[i:])
else:
int2 = 0
# Write new values to buffer variable
Bgrid_new[i] = C * (int1 + int2)
# Calculate relative error (take the maximum)
maxError = np.max(np.abs((Bgrid - Bgrid_new)/Bgrid_new))
# Copy values from buffer to working variable
Bgrid = np.copy(Bgrid_new)
# If the change in the last iteration was small enough, stop
if (maxError < 10**-4):
break
print "Finished in", j+1, "iterations, relative error:", maxError
return p2grid, Bgrid/np.sqrt(L2)
t0 = time.clock()
p2grid, Bgrid = solver()
print "Time:", time.clock() - t0, " seconds"
I am wondering if there a speedups possible:
- I am using
np.add.reduce
instead ofnp.sum
, which is about 0.2 secs faster on my system - I tried to move the code where
int1
andint2
are calculated to a separate function, no improvement here - Using
[... for enumerate(p2grid)]
instead of thefor
-loop seems not be faster either - The profiler says almost all time is still spent in the
reduce
method.
I can't see how this code can be optimized further (as most time is spent in a function of the NumPy lib), but I can't believe I wrote the most efficient code either.
If there are no optimizations possible at code level, is there anything else I should try or look into, like different interpreters or something like that?