Skip to content

Commit a1af8a5

Browse files
arm61arm61
authored
Merge pull request #21 from NHNahler/master
Introduction to scipy
2 parents 8b50967 + ce4b185 commit a1af8a5

File tree

1 file changed

+330
-0
lines changed

1 file changed

+330
-0
lines changed

‎notebooks/PiC_scipy.ipynb‎

Lines changed: 330 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,330 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"id": "9242e111",
6+
"metadata": {},
7+
"source": [
8+
"# Prerequisites\n",
9+
"- `matplotlib`\n",
10+
"- `numpy`\n",
11+
"- functions\n",
12+
"- f-strings"
13+
]
14+
},
15+
{
16+
"cell_type": "markdown",
17+
"id": "0dac4919",
18+
"metadata": {},
19+
"source": [
20+
"# Learning outcomes\n",
21+
"- Use physical and chemical constants from a library.\n",
22+
"- Convert between common units of temperature, pressure and energy.\n",
23+
"- Perform a linear regression and plot the result.\n",
24+
"- Perform a non-linear least-squares fit and plot the result."
25+
]
26+
},
27+
{
28+
"cell_type": "markdown",
29+
"id": "ad0fb146-f9ef-4d6c-9a04-ce3a7ad3ddd2",
30+
"metadata": {},
31+
"source": [
32+
"# Scipy\n",
33+
"A scientific Python library that among other functions provides:\n",
34+
"- scientific constants\n",
35+
"- advanced linear regression\n",
36+
"- non-linear least squares curve fitting"
37+
]
38+
},
39+
{
40+
"cell_type": "markdown",
41+
"id": "b33bf40b",
42+
"metadata": {},
43+
"source": [
44+
"## Use of scientific constants\n",
45+
"In any calculations or script you should define scientific constants with a clearly named variable rather than typing the number directly in your calculation. It is best practice to use the SI units for the constant."
46+
]
47+
},
48+
{
49+
"cell_type": "code",
50+
"execution_count": null,
51+
"id": "1731cb41",
52+
"metadata": {},
53+
"outputs": [],
54+
"source": [
55+
"# Ideal gas constant\n",
56+
"R = 8.314 # in J K^-1 mol^-1\n",
57+
"# The volume of 1 mol of ideal gas at a pressure of 1 bar and a temperature of 298 K.\n",
58+
"p = 100_000 # in Pa\n",
59+
"n = 1 # in mol\n",
60+
"T = 298 # in K\n",
61+
"print(f\"The volume of 1 mol of ideal gas at standard pressure and 298 K is {1000*n*R*T/p:.2f} L.\")"
62+
]
63+
},
64+
{
65+
"cell_type": "markdown",
66+
"id": "3e6b77c3",
67+
"metadata": {},
68+
"source": [
69+
"## Scientific constants with `scipy.constants`\n",
70+
"Import the `constants` package from the `scipy` library as a whole or import the specific constant(s) that you need."
71+
]
72+
},
73+
{
74+
"cell_type": "code",
75+
"execution_count": null,
76+
"id": "ec42bf6b",
77+
"metadata": {},
78+
"outputs": [],
79+
"source": [
80+
"from scipy.constants import R,c\n",
81+
"print(f\"The ideal gas constant is R = {R} J K^-1 mol^-1)\")\n",
82+
"print(f\"The speed of light is c = {c} m s^-1\")"
83+
]
84+
},
85+
{
86+
"cell_type": "code",
87+
"execution_count": null,
88+
"id": "62208e55",
89+
"metadata": {},
90+
"outputs": [],
91+
"source": [
92+
"import scipy.constants as c\n",
93+
"print(f\"The Boltzmann constant is k_B = {c.k} J K^-1\")\n",
94+
"print(f\"Planck's constant is h = {c.h} J s\")"
95+
]
96+
},
97+
{
98+
"cell_type": "markdown",
99+
"id": "59a13bef",
100+
"metadata": {},
101+
"source": [
102+
"## Exercise\n",
103+
"Use the constants from the `scipy.constants` package to convert the energy of 1 kJ mol<sup>-1</sup> into\n",
104+
"- eV\n",
105+
"- cm<sup>-1</sup>"
106+
]
107+
},
108+
{
109+
"cell_type": "code",
110+
"execution_count": null,
111+
"id": "5b030b98",
112+
"metadata": {},
113+
"outputs": [],
114+
"source": []
115+
},
116+
{
117+
"cell_type": "markdown",
118+
"id": "451e15a4",
119+
"metadata": {},
120+
"source": [
121+
"You can find a whole list of physical constants included in the `scipy.constants` package at https://docs.scipy.org/doc/scipy/reference/constants.html"
122+
]
123+
},
124+
{
125+
"cell_type": "markdown",
126+
"id": "c51b2c98",
127+
"metadata": {},
128+
"source": [
129+
"## Advanced linear regression with `scipy.stats.linregress`\n",
130+
"The linear regression method provided by the `scipy` package provides some advanced statistical parameters, those that users may be familiar with from Microsoft Excel's LINEST:\n",
131+
"- Slope\n",
132+
"- Intercept\n",
133+
"- R value\n",
134+
"- P value\n",
135+
"- standard error in the slope\n",
136+
"- standard error in the intercept"
137+
]
138+
},
139+
{
140+
"cell_type": "code",
141+
"execution_count": null,
142+
"id": "36ba6ccf",
143+
"metadata": {},
144+
"outputs": [],
145+
"source": [
146+
"# Define an x,y dataset that follows a linear trend\n",
147+
"import numpy as np\n",
148+
"x = np.arange(0,10)\n",
149+
"y = (3.0 + np.random.rand()) * (x + np.random.rand(10))"
150+
]
151+
},
152+
{
153+
"cell_type": "code",
154+
"execution_count": null,
155+
"id": "5d715e20",
156+
"metadata": {},
157+
"outputs": [],
158+
"source": [
159+
"# Perform linear regression\n",
160+
"from scipy.stats import linregress\n",
161+
"result = linregress(x,y)\n",
162+
"print(f\"The linear regression parameters are: {result}\")"
163+
]
164+
},
165+
{
166+
"cell_type": "markdown",
167+
"id": "c0d7c80a",
168+
"metadata": {},
169+
"source": [
170+
"### You can use `matplotlib` to visualise the data and the linear regression\n",
171+
"Each of the elements of the fitting results can be accessed by appending the name, separated by a dot from the regression result e.g. `result.slope` for the slope or `result.intercept` for the intercept."
172+
]
173+
},
174+
{
175+
"cell_type": "code",
176+
"execution_count": null,
177+
"id": "5ad05330",
178+
"metadata": {},
179+
"outputs": [],
180+
"source": [
181+
"import matplotlib.pyplot as plt\n",
182+
"plt.scatter(x,y,label=\"Data\")\n",
183+
"plt.plot(x,x*result.slope + result.intercept,label=\"Regression\")\n",
184+
"plt.title(\"Linear regression\")\n",
185+
"plt.xlabel(\"x\")\n",
186+
"plt.ylabel(\"y\")\n",
187+
"plt.legend()\n",
188+
"plt.show()"
189+
]
190+
},
191+
{
192+
"cell_type": "markdown",
193+
"id": "4b1acf38",
194+
"metadata": {},
195+
"source": [
196+
"## Exercise"
197+
]
198+
},
199+
{
200+
"cell_type": "markdown",
201+
"id": "4df88078",
202+
"metadata": {},
203+
"source": [
204+
"## Non-linear least-squares curve fitting with `scipy.optimize.curve_fit`\n",
205+
"While as Chemists we like to linearise data to then apply a linear regression, sometimes data is in a format that cannot be easily linearised and requires to fit the non-linear data with the appropriate function."
206+
]
207+
},
208+
{
209+
"cell_type": "code",
210+
"execution_count": null,
211+
"id": "e9a5c00a",
212+
"metadata": {},
213+
"outputs": [],
214+
"source": [
215+
"# Import required packages\n",
216+
"from scipy.optimize import curve_fit\n",
217+
"import numpy as np\n",
218+
"import matplotlib.pyplot as plt\n",
219+
"from scipy.constants import pi,k,u # Pi, Boltzmann constant, atomic mass unit"
220+
]
221+
},
222+
{
223+
"cell_type": "markdown",
224+
"id": "4cbc0b5c",
225+
"metadata": {},
226+
"source": [
227+
"### Maxwell-Boltzmann distribution\n",
228+
"We will determine the mean speed and temperature of argon."
229+
]
230+
},
231+
{
232+
"cell_type": "code",
233+
"execution_count": null,
234+
"id": "a3ecac55",
235+
"metadata": {},
236+
"outputs": [],
237+
"source": [
238+
"# Generate some data with noise\n",
239+
"T = 300 # in K\n",
240+
"steps = 50\n",
241+
"v = np.linspace(1,1000,steps)\n",
242+
"m = 39.948*u # mass of argon in atomic mass units\n",
243+
"F = ((m/(2*pi*k*T)**1.5 * 4*pi*v**2 * np.exp(-m*v**2/(2*k*T)))/1e9+np.random.rand(steps))*1e9\n",
244+
"plt.scatter(v,F)"
245+
]
246+
},
247+
{
248+
"cell_type": "code",
249+
"execution_count": null,
250+
"id": "5cb2ddc1",
251+
"metadata": {},
252+
"outputs": [],
253+
"source": [
254+
"def func(v,T,m):\n",
255+
" \"\"\"Fitting function for a Maxwell-Boltzmann distribution.\"\"\"\n",
256+
" return m/(2*pi*k*T)**1.5 * 4*pi*v**2 * np.exp(-m*v**2/(2*k*T))\n",
257+
"\n"
258+
]
259+
},
260+
{
261+
"cell_type": "code",
262+
"execution_count": null,
263+
"id": "4030e65a",
264+
"metadata": {},
265+
"outputs": [],
266+
"source": [
267+
"# Performing the nonlinear fit which returns the fitting parameters \n",
268+
"# as a tuple (which is like an immutable list) and the covariance matrix.\n",
269+
"#\n",
270+
"# In the fit we assume the element also as an unknown i.e. as a fitting parameter together with the temperature.\n",
271+
"popt, pcov = curve_fit(func, v, F, p0=[250,1e-25])\n"
272+
]
273+
},
274+
{
275+
"cell_type": "code",
276+
"execution_count": null,
277+
"id": "33b4692d",
278+
"metadata": {},
279+
"outputs": [],
280+
"source": [
281+
"plt.scatter(v,F,label=\"data\")\n",
282+
"plt.plot(v,func(v,*popt),label=f\"fit with T = {popt[0]:.1f} K, m = {popt[1]/u:.1f} amu\")\n",
283+
"plt.title(\"Fitted Maxwell-Boltzmann distribution\")\n",
284+
"plt.xlabel(r\"speed / m s$^{-1}$\")\n",
285+
"plt.ylabel(\"Intensitiy / arb. units\")\n",
286+
"plt.legend()\n",
287+
"plt.show()"
288+
]
289+
},
290+
{
291+
"cell_type": "markdown",
292+
"id": "d39c119f",
293+
"metadata": {},
294+
"source": [
295+
"# TODO\n",
296+
"- Add example for unit conversion with `scipy.constants`\n",
297+
"- Add exercises to all sections (unit conversion, linear regression possibly based on first linearising data, non-linear fit)\n",
298+
"- add some questions / quizzes\n",
299+
"- add a few debugging examples where suitable"
300+
]
301+
},
302+
{
303+
"cell_type": "markdown",
304+
"id": "a73b7e31",
305+
"metadata": {},
306+
"source": []
307+
}
308+
],
309+
"metadata": {
310+
"kernelspec": {
311+
"display_name": ".venv",
312+
"language": "python",
313+
"name": "python3"
314+
},
315+
"language_info": {
316+
"codemirror_mode": {
317+
"name": "ipython",
318+
"version": 3
319+
},
320+
"file_extension": ".py",
321+
"mimetype": "text/x-python",
322+
"name": "python",
323+
"nbconvert_exporter": "python",
324+
"pygments_lexer": "ipython3",
325+
"version": "3.11.2"
326+
}
327+
},
328+
"nbformat": 4,
329+
"nbformat_minor": 5
330+
}

0 commit comments

Comments
 (0)