Updating Chapter 9

Chapter_9/capacitor.py

@@ -36,11 +36,11 @@ distance = 0.1
 
 potential = V(0,0,distance)
 theory = distance*8*math.pi
-print("In the center : potential = {} - theory = {};  diff = {:.2f}%".format(potential,theory, 200*(potential-theory)/(potential + theory)))
+print("In the center : potential = {} - theory = {};  diff = {:.2f}%".format(potential,theory, 100*(potential-theory)/(theory)))
 
 
-potential = V(10,10,distance)
-print("In the corner : potential = {} - theory = {};  diff = {:.2f}%".format(potential,theory, 200*(potential-theory)/(potential + theory)))
-
+potential = V(9.5,9.5,distance)
+print("In the corner : potential = {} - theory = {};  diff = {:.2f}%".format(potential,theory, 100*(potential-theory)/(theory)))
 
 
+print("The two values are different because the assumption that we have a finite uniformly charged conductor is not physical.")

Chapter_9/ex1.py

+import random
+import math
+
+
+class my_random_vector:
+    def __init__(self,N_dim,min_value = -5, max_value = 5):
+        self.vec = [min_value + (max_value-min_value)*random.random() for i in range(N_dim)]
+
+    def norm_squared(self):
+        nn = 0
+        for x in self.vec:
+            nn+= x**2
+        return nn
+
+def pdf(vec):
+    return math.exp(-0.5*vec.norm_squared())
+
+
+def integrate_normal(n_dim, n_tests, n_sigma = 5):
+    n_tests = int(n_tests)
+    integral = 0
+    for test in range(n_tests):
+        integral += pdf(my_random_vector(n_dim,-n_sigma,n_sigma))
+    return 1./n_tests * integral * (2*n_sigma)**n_dim
+
+
+
+for n_dim in [1,2,5,15]:
+    mc = integrate_normal(n_dim,1e5,3)
+    truth = math.sqrt((2*math.pi)**n_dim)
+    print("N_dim = {}".format(n_dim))
+    print("Truth = {}, Monte-Carlo = {}".format(truth,mc))
+    print("Difference : {:.2f}%".format(100*(truth-mc)/truth))

Chapter_9/integrator.py

@@ -32,7 +32,7 @@ def TrapezeIntegrator(function, x_min, x_max, precision=1e-9, step_min = 3, step
     return integral
 
 
-def RombergIntegrator(function, x_min, x_max, precision=1e-9, step_min = 3, step_max = 20, silent = True):
+def RombergIntegrator(function, x_min, x_max, precision=1e-9, step_min = 3, step_max = 25, silent = True):
     
     r = [TrapezeStep(function, x_min, x_max, 1, 0)]
     
@@ -45,8 +45,10 @@ def RombergIntegrator(function, x_min, x_max, precision=1e-9, step_min = 3, step
         if not silent:
             print("Step {} : precision of {}".format(k, error))
         r = new_r
-        if error < precision and k > step_min:
+        if error < precision and k >= step_min:
             break
+    if k == step_max-1:
+        print("Warning, integral didn't converge after {} iterations with a precision of {}.".format(k, error))
     return r[-1]
     #self.integrals.append(self.r[-1])
 
@@ -55,9 +57,54 @@ def RombergIntegrator(function, x_min, x_max, precision=1e-9, step_min = 3, step
 
 
 if __name__ == "__main__":
-    def my_func(x):
-        return math.exp(-0.5*x**2)
+    print("### Trying with a gaussian function. The integral is expected to be {}".format(math.sqrt(2*math.pi)))
+    print("gaus [-5,5] : {}".format(RombergIntegrator(lambda x: math.exp(-0.5*x**2) , -5, 5 )))
+    print("gaus [-100,100] : {}".format(RombergIntegrator(lambda x: math.exp(-0.5*x**2), -100, 100 )))
+    print("\n\n### Trying with sin^2(x), should converge towards pi={}".format(math.pi))
+    print("Sin^2 with no minimum steps : {}".format(RombergIntegrator(lambda x : math.sin(x)**2, -math.pi, math.pi ,step_min = 0)))
+    print("This first integral converged towards a wrong value.")
+    print("It was expected since the first steps of Romberg are :")
+    print(" * evaluate R11 = (f(b) + f(a)) * (b-a)/2 = 0 since sin²(-pi) = sin²(pi) = 0")
+    print(" * evaluate R21 = previous_integral/2 + f( (b-a)/2 )/2 = 0/2 + 0/2 = 0")
+    print(" * compute  R22 = 1/(4^2 -1) * (4^2 * R21 - R11) = 0")
+    print("Since the first and second steps yield the same integral (0), the difference is < 1e-9 and the process stops")
+    print("This same integral converges correctly if we ask at least 3 steps: ")
+    print("Sin^2 : {}".format(RombergIntegrator(lambda x : math.sin(x)**2, -math.pi, math.pi)))
+    
+    
+    def test_1(x):
+        if x == 0:
+            return 0
+        return 1./math.sqrt(abs(x))
+    
+    print("\n\n### 1/sqrt(|x|) should converge towards 2")
+    print("1/sqrt(|x|) : {}".format(RombergIntegrator(test_1, -1, 1)))
+    print("We can see that it is indeed starting to converge towards the correct value, but very slowly...")
+    
+    
+    def test_2(x):
+        if x == 0:
+            return 0
+        return 1./x**2 - 1
+    
+    print("\n\n### 1/x^2 should diverge")
+    print("1/x^2 - 1 : {}".format(RombergIntegrator(test_2, -1 , 1 )))
+    print("And it is indeed the case")
+    
 
-    ii = RombergIntegrator(my_func, -5, 5 , silent = False)
-    print (ii)
+    def regular(x):
+        r  = x - int(x)
+        if r == 0:
+            return 0
+        return r**-2   
+    
+    print("\n\n### Finally, with the weird 'regular' function :")
+    print("regular [-256; 256] : {}".format(RombergIntegrator(regular, -256, 256)))
+    print("It converges towards 0 if we integrate on the whole range")
+    print("regular [0,1] : {}".format(RombergIntegrator(regular, 0, 1)))
+    print("But diverges when computing over [0,1]")
+    print("This is expected since :")
+    print("The first ten steps of Romberg will always evaluate on integer numbers yielding the same '0' value, thus returning a null integral")
+    print("But we know that it diverges if the integral is computed correctly (on the [0,1] interval for example)")
+