refer scipy.integrate in ex.3 hint

src/4_sommerfeld_model.md

@@ -354,7 +354,7 @@ A hypothetical metal has a Fermi energy $\epsilon_F = 5.2 \, \mathrm{eV}$ and a
 4. Now, find this difference in energy by calculating the integral found in 1 numerically. Compare your result with 3.
 
 	??? hint
-		 You can write a python script to compute integral using midpoint rule.
+		 You can do numerical integration in python with [`scipy.integrate.quad(func, xmin, xmax)`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html)
 
 5. Calculate the heat capacity for $T = 1000 \, \mathrm{K}$ in eV/K.
 6. Numerically compute the heat capacity by approximating the derivative of energy difference found in 4 with respect to $T$. To this end, make use of the fact that $$\frac{dy}{dx}=\lim_{\Delta x \to 0} \frac{y(x + \Delta x) - y(x - \Delta x)}{2 \Delta x}.$$ Compare your result with 5.

src/4_sommerfeld_model_solutions.md

@@ -101,6 +101,7 @@ kB = 8.617343e-5
 T = 1000 #kelvin
 
 import numpy as np
+from scipy import integrate
 
 np.seterr(over='ignore')
 
@@ -116,13 +117,12 @@ def g(E):
 def integral(E, T):
     return f(E, T)*g(E)*E
 
-## Solve integral using midpoint rule
-elements = np.linspace(0, 1e6, 1e7) #Free to choose the number of discrete steps
-dE = np.sum(integral(elements, T))*((1e6 - 0)/len(elements)) - 0.8e10 * 5.2**(5./2)
+## Solve integral numerically using scipy's integrate
+dE = integrate.quad(integral, 0, 1e1, args=(T))[0] - 0.8e10 * 5.2**(5./2)
 
 dT = 0.001
-dEplus = np.sum(integral(elements, 1000+dT))*((1e6 - 0)/len(elements)) - 0.8e10 * 5.2**(5./2)
-dEmin = np.sum(integral(elements, 1000-dT))*((1e6 - 0)/len(elements)) - 0.8e10 * 5.2**(5./2)
+dEplus = integrate.quad(integral, 0, 1e1, args=(T+dT))[0] - 0.8e10 * 5.2**(5./2)
+dEmin = integrate.quad(integral, 0, 1e1, args=(T-dT))[0] - 0.8e10 * 5.2**(5./2)
 
 CV = (dEplus - dEmin) / (2*dT);