The problem deals with the calculation of the deflection angle θ of two colliding structureless bodies under the effect of the potential V(r) (with r being the distance between the two bodies) without integrating the related Hamilton equations.
Such solution is obtained by transforming the 12 dimensional two body problem into that of a single particle of reduced mass µ using the center of mass frame and imposing the conservation of Energy E and the total angular momentum (see A. Laganà and G.A. Parker, Basic theory and computing for understanding chemical reactions, Springer (in the press)).
The expression for θ is:
in which b is the impact parameter.
To follow, we compare first the analytical solution:
for the potential $V(r) = \frac{1}{r}$ with the numerical quadrature of the integral given above.
Next we can see the corresponding Python code whose dependence on the impact parameter:
In [1]:
%matplotlib inline
from matplotlib import pyplot as plt
from matplotlib import patches as mpatches
from math import sqrt, degrees, asin
MAX_DIVISIONS = 21
MAX_ITERATIONS = 21
PRECISION_DELTA = 0.00000002
DELTA = 16.0
def analytic_evaluation(var_b, var_e):
"""Analytic function of theta integral."""
return 2.0 * asin(1 / (sqrt(1 + 4 * var_b ** 2 * var_e ** 2)))
def integral_func(var_r, var_b, var_e):
"""Integral to calculate."""
potential = lambda var_r: 1.0 / var_r
try:
return 1.0 / (
var_r**2 * sqrt(1 - (var_b**2 / var_r**2) - (
potential(var_r) / var_e)))
except ValueError:
return 0.0
except ZeroDivisionError:
return 0.0
def finite_integral(var_up, var_to, var_b, var_e):
"""Finite Riemann integral."""
prev_sum = 0.0
for i in range(MAX_DIVISIONS):
step = (var_to-var_up) / 2**i
partial_sum = 0.0
cur_step = var_up + step
while cur_step <= var_to:
f_x = step * integral_func(cur_step - (step/2.0), var_b, var_e)
partial_sum += f_x
cur_step += step
if prev_sum != 0.0 and abs(prev_sum - partial_sum) < PRECISION_DELTA:
break
prev_sum = partial_sum
return partial_sum
def integral_to_infinite(var_up, var_b, var_e):
"""Riemann integral open at right plus final multiplication."""
prev_res = 0.0
for i in range(MAX_ITERATIONS):
var_to = var_up + DELTA
partial_res = finite_integral(var_up, var_to, var_b, var_e)
if prev_res != 0.0 and partial_res != 0.0\
and partial_res < PRECISION_DELTA:
break
prev_res += partial_res
var_up = var_to
return (2.0 * var_b) * prev_res
var_a = 0.0
var_e = 0.1
points_x = []
points_y = []
points_y_a = []
print("Calculus may take some time...")
for var_b in range(0, 16):
points_x.append(var_b)
final_theta = 180.0 - degrees(abs(integral_to_infinite(var_a, var_b, var_e)))
analytic = degrees(analytic_evaluation(var_b, var_e))
print("b = {}\ttheta = {}\tanalytic = {}".format(var_b, final_theta, analytic))
points_y.append(final_theta)
points_y_a.append(analytic)
plt.ylabel('teta')
plt.xlabel('b')
plt.plot(points_x, points_y, 'r--')
plt.plot(points_x, points_y, 'ro')
plt.plot(points_x, points_y_a, 'b--')
plt.plot(points_x, points_y_a, 'bo')
red_line = mpatches.Patch(color='red', label='theta value')
blue_line = mpatches.Patch(color='blue', label='analytic value')
plt.legend(handles=[red_line, blue_line])
plt.grid(True)
plt.figure()
Out[1]:
The implementation of the program gave us values similar to the analytic function.
The major problems arouse from:
For the generic potential V(r) not having an analytical solution we did not consider to use the above described Python script because its performance is low due to the limitation of Python VM and of its interpreted nature.
We considered instead to use a combination of the C language and Python. This is also the first step to implement the code using C MPI API for parallel runs.
First of all the C code:
#ifdef _WIN32
#define _USE_MATH_DEFINES
#endif
#include <stdio.h>
#include <math.h>
const size_t MAX_DIVISIONS = 21;
const size_t MAX_ITERATIONS = 21;
const double PRECISION_DELTA = 0.00000001;
const double DELTA = 16.0;
double potential(double r)
{
return 1.0/r;
}
double function(double r, double b, double E)
{
return 1.0 / (pow(r, 2) * sqrt(1 - (pow(b, 2) / pow(r, 2)) - (potential(r) / E)));
}
double finite_integral(double up, double to, double b, double E)
{
double cur_step,
step,
prev_sum,
partial_sum,
f_x;
size_t i;
prev_sum = 0.0;
for(i = 0; i != MAX_DIVISIONS; ++i)
{
step = (to-up)/pow(2, i);
partial_sum = 0.0;
for(cur_step = up + step; cur_step <= to; cur_step += step) {
f_x = step * function(cur_step - (step/2.0), b, E);
// when f_x is nan
// if f_x is nan, f_x != f_x will be true
if(f_x == f_x) {
partial_sum = partial_sum + f_x;
}
}
if(prev_sum != 0.0 && fabs(prev_sum - partial_sum) < PRECISION_DELTA) break;
prev_sum = partial_sum;
}
return partial_sum;
}
double integral_to_infinite(double a, double b, double E)
{
double to,
prev_res,
partial_res;
size_t i;
prev_res = 0.0;
for(i = 1; i != pow(2, MAX_ITERATIONS); ++i)
{
to = a + DELTA;
partial_res = finite_integral(a, to, b, E);
if(prev_res != 0.0 && partial_res != 0.0 && partial_res < PRECISION_DELTA)
{
prev_res += partial_res;
break;
}
prev_res += partial_res;
a = to;
}
return (2.0 * b) * prev_res;
}
double to_degrees(double radians) {
return 180.0 - radians * (180.0 / M_PI);
}
The approach is the same of the Python code given above. Now we can see how to implement this code with Python:
from cffi import FFI
from matplotlib import pyplot as plt
from matplotlib import patches as mpatches
from ctypes import c_double
from math import asin, sqrt, degrees
def analytic_evaluation(var_b, var_e):
"""Analytic function of theta integral."""
return 2.0 * asin(1 / (sqrt(1 + 4 * float(var_b) ** 2 * float(var_e) ** 2)))
def range_f(start, end, step):
"""Create a list of double values."""
end = c_double(end)
ffi = FFI()
start = c_double(start)
while start.value <= end.value:
yield ffi.cast("double", start.value)
start.value = start.value + c_double(step).value
def main():
"""Program main."""
ffi = FFI()
ffi.cdef("""
double integral_to_infinite(double a, double b, double E);
double to_degrees(double radians);
""")
lib = ffi.verify("""
#include "tetaQuad.h"
""", include_dirs=["."], extra_compile_args=["-O3"])
var_e = ffi.cast("double", 0.1)
var_a = ffi.cast("double", 0.0)
points_x = []
points_y = []
points_y_a = []
print("Calculus may take some time...")
for var_b in range_f(0.0, 100, 0.5):
points_x.append(var_b)
rad = lib.integral_to_infinite(var_a, var_b, var_e)
theta = lib.to_degrees(abs(rad))
analytic = degrees(analytic_evaluation(var_b, var_e))
print("b = {}\ttheta = {}\tanalytic = {}".format(
float(var_b), theta, analytic))
points_y.append(theta)
points_y_a.append(analytic)
plt.ylabel('teta')
plt.xlabel('b')
plt.plot(points_x, points_y, 'ro')
plt.plot(points_x, points_y_a, 'b--')
red_line = mpatches.Patch(color='red', label='theta value')
blue_line = mpatches.Patch(color='blue', label='analytic value')
plt.legend(handles=[red_line, blue_line])
plt.grid(True)
plt.show()
if __name__ == '__main__':
main()
As we can see now the performance that we reached with the C implementation gave us a better chart in less than a minute. Obviously this is a start point to improve the code with a parallel implementation, searching for the tasks that can run concurrently.
Here we consider the Lennard-Jones potential (attractive-repulsive):
The chart values are almost as we expected, and still there is room for improvement. For example, we can improve precision by refining the quadrature method especially when using more complex kinds of potential functions.
The unit of measure used are submultiples of the m and s, respectively nanometer ($10^{-9}m$) and millisecond ($10^{-3}s$). For the potential we can misure it in millielectron ($10^{-3}eV, meV$).