Adding integration

Chapter_8/exam_simul.py

@@ -26,13 +26,8 @@ n_worse_sr = 0
 n_tries = 5000
 for i in range(n_tries):
     ss = simul_one_year(results_career,30)
-    av = ss.get_average()
-    sr = get_success_rate(ss)
-
-    if av < one_year_av:
-        n_worse_av +=1
-    if sr < one_year_sr:
-        n_worse_sr +=1
+    n_worse_av += ss.get_average() < one_year_av
+    n_worse_sr += get_success_rate(ss) < one_year_sr
 
 print("Worse average in {} % of simulations.".format(100*n_worse_av*1./n_tries))
 print("Worse success rate in {} % of simulations.".format(100*n_worse_sr*1./n_tries))

Chapter_8/tp_lphys1201.py

@@ -143,6 +143,7 @@ def sweep_alpha(n_students, wait_thresh, n_simul = 1000, n_teachers = 2, av_stud
 
 
 simulate_many(10000 , 13 , 2 , 40 )
+simulate_many(10000 , 13 , 1 , 40 )
 
 simulate_many(10000 , 26 , 2 , 40 )
 

Chapter_9/capacitor.py

+from integrator import RombergIntegrator as integral
+import math
+
+
+def d_field(rx,ry,rz,Q):
+    den = rx**2 + ry**2 + rz**2
+    if den == 0:
+        return 0
+
+    return Q * rz / math.pow(den, 1.5)
+
+
+def I_x(x, ry, rz, Q, L):
+
+    def ix(rx):
+        return d_field(rx,ry,rz,Q)
+
+    return integral(ix, x-L, x+L, precision = 1e-2)
+
+def E_z(x,y,rz, Q = 1, L = 10):
+
+    def iy(ry):
+        return I_x(x,ry,rz,Q,L)
+
+    return integral(iy,y-L, y+L, precision = 1e-2)
+
+
+def V(x,y, z0 = 1, Q = 1, L = 10):
+    def vz(rz):
+        return 2*E_z(x,y,rz,Q, L)
+    return integral(vz,0,2*z0,precision = 1e-2,silent = False)
+
+
+distance = 0.2
+
+
+potential = V(0,0,distance/2)
+theory = distance*4*math.pi
+print("In the center : potential = {} - theory = {};  diff = {:.2f}%".format(potential,theory, 200*(potential-theory)/(potential + theory)))
+
+
+potential = V(10,10,distance)
+print("In the corner : potential = {} - theory = {};  diff = {:.2f}%".format(potential,theory, 200*(potential-theory)/(potential + theory)))
+
+
+

Chapter_9/integrator.py

+import math
+        
+def TrapezeStep(function, x_min, x_max, k, old_int):
+    if k == 1:
+        return (function(x_min)+function(x_max))*(x_max-x_min)/2
+    else:
+        n_steps = 2**(k-2)
+        h = (x_max - x_min)/n_steps
+        ii = 0
+        x = x_min + h/2
+        for i in range(n_steps):
+            ii += function(x)
+            x  += h
+        return old_int/2 + h*ii/2
+
+
+def TrapezeIntegrator(function, x_min, x_max, precision=1e-9, step_min = 3, step_max = 25, silent = True):
+    
+    previous_integral = 0
+    integral = 0
+    error = precision+1
+    
+    for step in range(1,step_max):
+        integral = TrapezeStep(function, x_min, x_max, step, previous_integral)
+        error = abs(previous_integral - integral)
+        previous_integral = integral
+        if not silent:
+            print("Step {} : precision of {}".format(step, error))
+        if error < precision and step > step_min:
+            break
+    
+    return integral
+
+
+def RombergIntegrator(function, x_min, x_max, precision=1e-9, step_min = 3, step_max = 20, silent = True):
+    
+    r = [TrapezeStep(function, x_min, x_max, 1, 0)]
+    
+    for k in range(2,step_max):
+        new_r = [ TrapezeStep(function, x_min, x_max, k, r[0]) ]
+        for m in range(1,k):
+            new_r.append(1/(4**m - 1) * (4**m * new_r[m-1] - r[m-1]))
+            
+        error = abs(r[-1] - new_r[-1])
+        if not silent:
+            print("Step {} : precision of {}".format(k, error))
+        r = new_r
+        if error < precision and k > step_min:
+            break
+    return r[-1]
+    #self.integrals.append(self.r[-1])
+
+
+
+
+
+if __name__ == "__main__":
+    def my_func(x):
+        return math.exp(-0.5*x**2)
+
+    ii = RombergIntegrator(my_func, -5, 5 , silent = False)
+    print (ii)
+