From b5ff3688db0d656b49cb0c74fd41cc3714c38cc8 Mon Sep 17 00:00:00 2001
From: Martin Delcourt <martin.delcourt@uclouvain.be>
Date: Fri, 22 Nov 2019 09:19:12 +0100
Subject: [PATCH] Updating Chapter 9
---
Chapter_9/capacitor.py | 8 +++---
Chapter_9/ex1.py | 33 +++++++++++++++++++++++
Chapter_9/integrator.py | 59 ++++++++++++++++++++++++++++++++++++-----
3 files changed, 90 insertions(+), 10 deletions(-)
create mode 100644 Chapter_9/ex1.py
diff --git a/Chapter_9/capacitor.py b/Chapter_9/capacitor.py
index 8372cec..04c8ec3 100644
--- a/Chapter_9/capacitor.py
+++ b/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.")
diff --git a/Chapter_9/ex1.py b/Chapter_9/ex1.py
new file mode 100644
index 0000000..6b9d896
--- /dev/null
+++ b/Chapter_9/ex1.py
@@ -0,0 +1,33 @@
+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))
diff --git a/Chapter_9/integrator.py b/Chapter_9/integrator.py
index e7a568f..4534823 100644
--- a/Chapter_9/integrator.py
+++ b/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)")
+
--
GitLab