Average Force (Singularity-Free)

Applying Newton's law of gravitation, the magnitude of the force between two particles is:

$F(x_1,y_1,x_2,y_2) = \frac{1}{(x_2-x_1)^2 + (y_2-y_1)^2+1}$

The coordinates of each particle are uniformly distributed and the expected value of force can be computed as such:

$F_{av} = \int_{0}^{1} \int_{0}^{1} \int_{0}^{1} \int_{0}^{1} F(x_1,y_1,x_2,y_2) \ dx_1 \ dy_1 \ dx_2 \ dy_2 \approx 0.779256$

Integral evaluated using Wolfram-Alpha. Using a script of code, a step size of 0.01 does not converge to a reasonably accurate result and decreasing the step size would significantly increase the time taken for computation ( $\mathrm{O}(n^4)$ ) .

The proper way to solve this is with a quadruple integral, as shown in the solution by @Karan Chatrath . The downside, at least when solving numerically, is that a quadruple integral is very computationally expensive. To get a decent result, one would likely need between $10^8$ and $10^{12}$ loop iterations. By contrast, the Monte Carlo simulation method is very cheap. With just $10^6$ random trials, I can get a very good result. The Monte Carlo method requires the generation of random numbers with a uniform probability distribution. See code below:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

import math
import random

N = 10**6

Fav = 0.0

G = 1.0
m = 1.0

#########################################################

for j in range(0,N):

    x1 = random.uniform(0.0,1.0)
    y1 = random.uniform(0.0,1.0)
    z1 = 0.0

    x2 = random.uniform(0.0,1.0)
    y2 = random.uniform(0.0,1.0)
    z2 = 1.0

    Dx = x1 - x2
    Dy = y1 - y2
    Dz = z1 - z2

    Dsq = Dx**2.0 + Dy**2.0 + Dz**2.0

    F = G*m*m/Dsq

    Fav = Fav + F

#########################################################

Fav = Fav/N

print N
print Fav


#>>> 
#1000000
#0.779365934848
#>>> ================================ RESTART ================================
#>>> 
#1000000
#0.779299868788
#>>> ================================ RESTART ================================
#>>> 
#1000000
#0.779052645254
#>>> ================================ RESTART ================================
#>>> 
#1000000
#0.779407268843
#>>>