From 67838e3504abd8abbed38aff5312051be9f9a909 Mon Sep 17 00:00:00 2001
From: Sathish Kumar RK <rksathish09@gmail.com>
Date: Mon, 2 Mar 2020 15:06:07 +0000
Subject: [PATCH] refer scipy.integrate in ex.3 hint
---
src/4_sommerfeld_model.md | 2 +-
src/4_sommerfeld_model_solutions.md | 10 +++++-----
2 files changed, 6 insertions(+), 6 deletions(-)
diff --git a/src/4_sommerfeld_model.md b/src/4_sommerfeld_model.md
index 86c2390e..34096020 100644
--- a/src/4_sommerfeld_model.md
+++ b/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.
diff --git a/src/4_sommerfeld_model_solutions.md b/src/4_sommerfeld_model_solutions.md
index a520cc44..4c637fe5 100644
--- a/src/4_sommerfeld_model_solutions.md
+++ b/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);
--
GitLab