Mathematical Computing¶
The Nature of Mathematical Modeling (draft)¶
$ \def\%#1{\mathbf{#1}} \def\mat#1{\mathbf{#1}} \def\*#1{\vec{#1}} \def\ve#1{\vec{#1}} \def\ds#1#2{\cfrac{d #1}{d #2}} \def\dd#1#2{\cfrac{d^2 #1}{d #2^2}} \def\ps#1#2{\cfrac{\partial #1}{\partial #2}} \def\pp#1#2{\cfrac{\p^2 #1}{\p #2^2}} \def\p{\partial} \def\ba{\begin{eqnarray}} \def\ea{\end{eqnarray}} \def\eps{\epsilon} \def\del{\nabla} $
Open- vs closed-source¶
- Will use the former, for accessibility, portability, and extensibility
Languages¶
Python¶
- Combines best features of C++'s objects with Lisp's list processing and APL's data structure manipulation.
- Interpreted (interactive but slow), but can be accelerated.
- tutorial
Package management¶
- Virtual environments collect dependencies and minimize conflicts
- python -m venv environment_name
- source environment_name/bin/activate
- Proliferation of options
import this
Rust¶
- Fixes recurring language security and stability issues: buffer overflows, memory leaks, race conditions, ...
- Compiled, can target multiple platforms
- tutorial
Julia¶
- Developed as a high-level language with low-level performance
LLVM¶
- compiler infrastructure used by many languages
Web¶
Javascript/ECMAScript¶
- Browser scripting language with just-in-time (JIT) compilation
WebAssembly (WASM)¶
WebGPU¶
- Graphics and compute acceleration
Coding¶
Text editors¶
Integrated Development Environments (IDEs)¶
Large Language Models¶
- vibe coding
- describe rather than write code
- distill and apply documentation
- can make mistakes, infringe copyrights, ...
- Cursor, Codex, Copilot, Claude, Gemini, ...
- MCP
Notebooks¶
- Combine coding, development, documentation, distribution, collaboration
Jupyter¶
- cells
- code
- Markdown
- raw
- kernels
- execute code in cells
- JupyterLab
- JupyterHub
- user group server
- The Littlest JupyterHub
- Binder
- computing environment server
- Pretty Jupyter
- report formatting
- metadata
- Book
- book formatting
- Lite
- runs entierly in browser, using WebAssembly
- CoLab
- freemium cloud service
- vscode-jupyter
- VS code extension
- jupyterlab-vim
- Vi key bindings
- magic commands
- %pip
- %run
- %%javascript
- ipywidgets
- user interface elements
import ipywidgets as widgets
def button_update():
global count
output.clear_output()
count += 1
print(f"button pressed {count} times")
def button_handler(self):
with output:
button_update()
button = widgets.Button(description='click me')
output = widgets.Output()
count = 0
button.on_click(button_handler)
display(button,output)
Button(description='click me', style=ButtonStyle())
Output()
Version Control¶
Git¶
- distributed development
- clone, pull, commit, push, branch, merge
- tutorial
- jupyterlab-git
Computer¶
Architectures¶
- CPU
- Small numbers (tens) of general-purpose cores
- GPU
- Large number (thousands) of special-purpose cores
- MCU
- Low-cost, low-power, low-pin-count processors for embedded computing
Systems¶
clusters¶
- Slurm
- cluster management and job scheduling
- Docker Swarm
- deploying containerized applications
- Kubernetes (K8s)
- scaling containerized applications
- Zero to JupyterHub with Kubernetes
clouds¶
Threads¶
- Lightweight processes that can be run independently, and can distribute work over cores
import threading,time
def task(delay):
print(f" second thread: delay {delay}s")
time.sleep(delay)
print(" second thread: hello")
thread = threading.Thread(target=task,args=(1,))
print("main thread: hello")
thread.start()
thread.join()
print("main thread: second thread finished")
main thread: hello second thread: delay 1s second thread: hello main thread: second thread finished
Benchmarking¶
pi calculation¶
- simple scatter-gather test $$ \ba \pi &=& 4\tan^{-1}(1) \\ &\approx& \sum_{i=1}^N \frac{0.5}{(i-0.75)(i-0.25)} \ea $$
Python¶
import time
NPTS = 10000000
print("\nPython pi calculation...\n")
pi = 0
start_time = time.time()
for i in range(1,(NPTS+1)):
pi += 0.5/((i-0.75)*(i-0.25))
end_time = time.time()
mflops = NPTS*5.0/(1.0e6*(end_time-start_time))
print("NPTS = %d, pi = %f"%(NPTS,pi))
print("time = %f, estimated MFlops = %f"%(end_time-start_time,mflops))
Python pi calculation... NPTS = 10000000, pi = 3.141593 time = 0.514727, estimated MFlops = 97.138940
C¶
$ gcc pi.c -o pi -lm -O3
$ ./pi
NPTS = 1000000000, pi = 3.141593
time = 0.686136, estimated MFlops = 7287.190035
$ gcc pi.c -o pi -lm -O3 -ffast-math
$ ./pi
NPTS = 1000000000, pi = 3.141593
time = 0.353395, estimated MFlops = 14148.484225
$ g++ threadpi.cpp -o threadpi -O3 -ffast-math -pthread
$ ./threadpi
npts: 100000000 nthreads: 16 pi: 3.14159
time: 0.085009 estimated MFlops: 94107.7
Packages¶
- High-level routines for mathematical computing
Numba¶
- JIT compiler for Python
from numba import njit
import time
NPTS = 10000000
print("\nNumba pi calculation ...\n")
import time
NPTS = 100000000
@njit(fastmath=True)
def calc():
pi = 0
for i in range(1,(NPTS+1)):
pi += 0.5/((i-0.75)*(i-0.25))
return pi
print("\nNumba serial version:")
pi = calc() # first call to compile the function
start_time = time.time()
pi = calc() # second call uses the cached compilation
end_time = time.time()
mflops = NPTS*5.0/(1.0e6*(end_time-start_time))
print("NPTS = %d, pi = %f"%(NPTS,pi))
print("time = %f, estimated MFlops = %f"%(end_time-start_time,mflops))
#
print("\nNumba parallel version")
from numba import prange
NPTS = 10000000000
@njit(parallel=True,fastmath=True)
def calc():
pi = 0
for i in prange(1,(NPTS+1)):
pi += 0.5/((i-0.75)*(i-0.25))
return pi
pi = calc() # first call to compile the function
start_time = time.time()
pi = calc() # second call uses the cached compilation
end_time = time.time()
mflops = NPTS*5.0/(1.0e6*(end_time-start_time))
print("NPTS = %d, pi = %f"%(NPTS,pi))
print("time = %f, estimated MFlops = %f"%(end_time-start_time,mflops))
Numba pi calculation ... Numba serial version: NPTS = 100000000, pi = 3.141593 time = 0.066294, estimated MFlops = 7542.165815 Numba parallel version NPTS = 10000000000, pi = 3.141593 time = 0.576868, estimated MFlops = 86674.968455
NumPy¶
- Mathematical computing routines for Python
import numpy as np
import time
NPTS = 10000000
print("\nNumPy pi calculation ...\n")
start_time = time.time()
i = np.arange(1,(NPTS+1),dtype=float)
pi = np.sum(0.5/((i-0.75)*(i-0.25)))
end_time = time.time()
mflops = NPTS*5.0/(1.0e6*(end_time-start_time))
print("NPTS = %d, pi = %f"%(NPTS,pi))
print("time = %f, estimated MFlops = %f"%(end_time-start_time,mflops))
NumPy pi calculation ... NPTS = 10000000, pi = 3.141593 time = 0.060566, estimated MFlops = 825.543138
Jax¶
- Python just-in-time compilation, parallel acceleration, automatic differentiation
print('Jax pi calculation ...\n')
#
import os,time
NCPU = os.cpu_count()
os.environ["XLA_FLAGS"] = f"--xla_force_host_platform_device_count={NCPU}"
import jax
jax.config.update('jax_platform_name','cpu')
import jax.numpy as jnp
# values for calculating pi
a = 0.5
b = 0.75
c = 0.25
# alternate compilation values to prevent caching
a0 = 0.6
b0 = 0.7
c0 = 0.2
#
print("compile Jax version:")
NPTS = 100000000
@jax.jit
def calcpi_jit(a,b,c):
i = jnp.arange(1,(NPTS+1),dtype=float)
pi = jnp.sum(a/((i-b)*(i-c)))
return pi
start_time = time.time()
pi = calcpi_jit(a0,b0,c0).block_until_ready()
end_time = time.time()
print("time = %f"%(end_time-start_time))
#
print("run Jax version:")
start_time = time.time()
pi = calcpi_jit(a,b,c).block_until_ready()
end_time = time.time()
mflops = NPTS*5.0/(1.0e6*(end_time-start_time))
print("NPTS = %d, pi = %f"%(NPTS,pi))
print("time = %f, estimated MFlops = %f"%(end_time-start_time,mflops))
#
print("\ncompile Jax parallel version:")
#
from jax.sharding import PartitionSpec as P,NamedSharding,AxisType
mesh = jax.make_mesh((NCPU,),('x'),axis_types=(AxisType.Explicit))
jax.set_mesh(mesh)
shard = jax.NamedSharding(mesh,jax.P('x'))
#
NPTS = 100000000*NCPU
@jax.jit
def calcpi_shard(a,b,c):
i = jnp.arange(1,(NPTS+1),dtype=float,device=shard)
pi = jnp.sum(a/((i-b)*(i-c)))
return pi
start_time = time.time()
pi = calcpi_shard(a0,b0,c0).block_until_ready()
end_time = time.time()
print("time = %f"%(end_time-start_time))
#
print(f"run Jax parallel version using {NCPU} CPUs:")
start_time = time.time()
pi = calcpi_shard(a,b,c).block_until_ready()
end_time = time.time()
mflops = NPTS*5.0/(1.0e6*(end_time-start_time))
print("NPTS = %d, pi = %f"%(NPTS,pi))
print("time = %f, estimated MFlops = %f"%(end_time-start_time,mflops))
Jax pi calculation ... compile Jax version: time = 0.130019 run Jax version: NPTS = 100000000, pi = 3.141592 time = 0.048594, estimated MFlops = 10289.286082 compile Jax parallel version: time = 0.366437 run Jax parallel version using 16 CPUs: NPTS = 1600000000, pi = 3.141592 time = 0.399800, estimated MFlops = 20010.001890
Jax pi calculation on H200 GPU ...
NPTS = 2000000000, pi = 3.141592
time = 0.001954, estimated MFlops = 5118750.305101
PyTorch, TensorFlow¶
- Widely-used machine learning packages
Scikit-learn¶
- High-level machine learning routines
Pandas¶
- Python data manipulation tools
OpenCV¶
- Standard package for computer vision
SciPy¶
- Scientific computing routines for Python
Sympy¶
- Symbolic computing for Python
- simplification, limits
import numpy as np
import sympy as sp
from IPython.display import display,Math
print('sqrt(8):')
print(f" NumPy: {np.sqrt(8)}")
print(f" SymPy: {sp.sqrt(8)}")
from sympy.abc import x
f = (x+1)
display(Math(f"f(x)={f}"))
display(Math(f"f(x)^{{10}} = {sp.latex(sp.expand(f**10))}"))
g = sp.lambdify(x,f**10,'numpy')
display(Math(f"f(1)^{{10}} = {g(1)}"))
sqrt(8): NumPy: 2.8284271247461903 SymPy: 2*sqrt(2)$\displaystyle f(x)=x + 1$ $\displaystyle f(x)^{10} = x^{10} + 10 x^{9} + 45 x^{8} + 120 x^{7} + 210 x^{6} + 252 x^{5} + 210 x^{4} + 120 x^{3} + 45 x^{2} + 10 x + 1$ $\displaystyle f(1)^{10} = 1024$
Libraries¶
- Low-level routines for mathematical computing
Numerical Reciptes¶
- Well-documented broad range of routines
BLAS (Basic Linear Algebra Subprograms)¶
- Widely-used core vector and matrix routines
MPI¶
- Widely-used distributed computing routines
CUDA¶
- NVIDIA GPU accelerated computing library
OpenCL¶
- Open library for accelerated computing
SYCL¶
- Open library for heterogeneous computing
Visualization¶
ipycanvas¶
- Notebook canvas drawing elements
import numpy as np
import time
from ipycanvas import Canvas,hold_canvas
import ipywidgets as widgets
from IPython.display import display
NPTS = 500
l = 15
x = np.linspace(-l,l,NPTS)
y = np.linspace(-l,l,NPTS)
x,y = np.meshgrid(x,y)
rcanvas = np.sqrt(x*x+y*y)
canvas = Canvas(width=NPTS,height=NPTS)
def canvas_button_handler(self):
for i in range(300):
with hold_canvas():
k = 2*i/300
z = 255*(1+np.cos(k*rcanvas))/2
image = np.dstack((z,z,z))
canvas.put_image_data(image,0,0)
time.sleep(0.005)
for i in range(300,0,-1):
with hold_canvas():
k = 2*i/300
z = 255*(1+np.cos(k*rcanvas))/2
image = np.dstack((z,z,z))
canvas.put_image_data(image,0,0)
time.sleep(0.005)
canvas_button = widgets.Button(description='animate')
canvas_button.on_click(canvas_button_handler)
display(canvas_button)
display(canvas)
Button(description='animate', style=ButtonStyle())
Canvas(width=500)
Matplotlib¶
- Wide range of data visualization routines
- Good control over graph layout
import matplotlib.pyplot as plt
import numpy as np
l = 15.5
x = np.arange(-l,l,0.2)
y = np.sin(x)/x
plt.plot(x,y)
plt.show()
#
# enable interactive features
#
%matplotlib ipympl
#
# imports
#
import ipywidgets as widgets
import matplotlib.pyplot as plt
import numpy as np
#
# set up interactive plot
#
l = 15.5
k = 1
xion = np.arange(-l,l,0.2)
yion = np.sin(k*xion)/(k*xion)
plt.ion()
fig,ax = plt.subplots()
fig.canvas.header_visible = False
line, = ax.plot(xion,yion)
plt.title('sin(kx)/(kx)')
plt.show()
#
# handle slider changes
#
def slider_handler(change):
k = change['new']
yion = np.sin(k*xion)/(k*xion)
line.set_ydata(yion)
fig.canvas.draw_idle()
#
# add slider
#
slider = widgets.FloatSlider(value=1,min=0.01,max=10.0,step=0.1,
description='adjust k:',
continuous_update=True,
orientation='horizontal',
readout_format='.1f',)
slider.observe(slider_handler,names='value')
display(slider)
FloatSlider(value=1.0, description='adjust k:', max=10.0, min=0.01, readout_format='.1f')
Plotly¶
- Wide range of data visualization routines
- Good interactivity
- Plotly Express
- high-level visualization interface
- Dash
- low-code application framework
import plotly.express as px
import numpy as np
l = 15.5
x = np.arange(-l,l,0.2)
y = np.sin(x)/x
pxfig = px.line(x=x,y=y)
pxfig.show()
import numpy as np
import plotly.graph_objects as go
import plotly.io as pio
pio.renderers.default = 'jupyterlab'
l = 15.5
gox = np.arange(-l,l,0.2)
goy = np.arange(-l,l,0.2)
(gox,goy) = np.meshgrid(gox,goy)
gor = np.sqrt(gox**2+goy**2)
goz = np.sin(gor)/gor
gofig = go.Figure(data=[go.Surface(z=goz,x=gox,y=goy,colorscale='solar')])
gofig.update_layout(width=750,height=750)
gofig.show()
import numpy as np
import plotly.graph_objects as go
import plotly.io as pio
import time
import ipywidgets as widgets
from IPython.display import display
pio.renderers.default = 'jupyterlab'
l = 15.5
gox = np.arange(-l,l,0.2)
goy = np.arange(-l,l,0.2)
(gox,goy) = np.meshgrid(gox,goy)
gor = np.sqrt(gox**2+goy**2)
goz = np.sin(gor)/gor
gowfig = go.FigureWidget()
gowfig.add_surface(z=goz,x=gox,y=goy,colorscale='solar')
gowfig.update_layout(width=750,height=750)
def go_button_handler(self):
k = 0
for i in range(50):
k += .05
goz = np.sin(k*gor)/(k*gor)
gowfig.update_traces(z=goz)
time.sleep(0.1)
gobutton = widgets.Button(description='animate')
gobutton.on_click(go_button_handler)
display(gobutton)
display(gowfig)
Button(description='animate', style=ButtonStyle())
FigureWidget({
'data': [{'colorscale': [[0.0, 'rgb(51, 19, 23)'], [0.09090909090909091,
'rgb(79, 28, 33)'], [0.18181818181818182, 'rgb(108,
36, 36)'], [0.2727272727272727, 'rgb(135, 47, 32)'],
[0.36363636363636365, 'rgb(157, 66, 25)'],
[0.45454545454545453, 'rgb(174, 88, 20)'],
[0.5454545454545454, 'rgb(188, 111, 19)'],
[0.6363636363636364, 'rgb(199, 137, 22)'],
[0.7272727272727273, 'rgb(209, 164, 32)'],
[0.8181818181818182, 'rgb(217, 192, 44)'],
[0.9090909090909091, 'rgb(222, 222, 59)'], [1.0,
'rgb(224, 253, 74)']],
'type': 'surface',
'uid': '858fae0d-0e54-4fca-95c7-d325c7202f0e',
'x': {'bdata': ('AAAAAAAAL8CamZmZmZkuwDQzMzMzMy' ... 'zMzMwtQPYyMzMzMy5AXJmZmZmZLkA='),
'dtype': 'f8',
'shape': '155, 155'},
'y': {'bdata': ('AAAAAAAAL8AAAAAAAAAvwAAAAAAAAC' ... 'mZmZkuQFyZmZmZmS5AXJmZmZmZLkA='),
'dtype': 'f8',
'shape': '155, 155'},
'z': {'bdata': ('G5hiHo5zaj+aIjhQlMSDP6qrTYLLTp' ... 't9m22cPymIY3L8nJY/icKxICtkkD8='),
'dtype': 'f8',
'shape': '155, 155'}}],
'layout': {'height': 750, 'template': '...', 'width': 750}
})
Three.js¶
- 3D graphics for JavaScript
- pythreejs
- Jupyter/Python bridge for ThreeJS
from pythreejs import *
import numpy as np
import time,math
import ipywidgets as widgets
from IPython.display import display
l = 15.5
nx = 100
x = np.linspace(-l,l,nx)
y = np.linspace(-l,l,nx)
x,y = np.meshgrid(x,y)
r = np.sqrt(x**2+y**2)
z = np.sin(r)/r
surf = SurfaceGeometry(z=list(z.flat),width=1,height=1,
width_segments=nx-1,height_segments=nx-1)
matl = MeshLambertMaterial(color='gold',side='DoubleSide')
mesh = Mesh(surf,matl)
camera = PerspectiveCamera(position=[0,-2,1],up=[0,1,0],
children=[DirectionalLight(color='white',position=[-10,0,0],intensity=1)])
scene = Scene(children=[mesh,camera,AmbientLight(color='white',intensity=0.1)])
renderer = Renderer(camera=camera,scene=scene,antialias=True,
controls=[OrbitControls(controlling=camera)],
width=500,height=500)
def three_button_handler(self):
for i in range(1,100):
k = 2*i/100
surf.z = list((np.sin(k*r)/(k*r)).flat)
time.sleep(.025)
three_button = widgets.Button(description='animate')
three_button.on_click(three_button_handler)
display(three_button)
display(renderer)
Button(description='animate', style=ButtonStyle())
Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', position=(-10.0, 0.0, 0.0), quater…
Sonification¶
- listening to data
- our ears are good at detecting features and patterns in sounds ### ipytone ### display
import numpy as np
from IPython.display import Audio
rate = 44100
duration = 5
fmin = 100
fmax = 500
t = np.linspace(0,duration,rate*duration)
data = np.sin(2*np.pi*(fmin+(fmax-fmin)*(t/duration))*t)
Audio(data,rate=rate)
References¶
- [Press:2007] Press, William H. Numerical recipes 3rd edition: The art of scientific computing. Cambridge university press, 2007.
- Very readable and very comprehensive survey of scientific computing
Problems¶
- Simulate/plot/render/... bouncing balls in a Jupyter notebook
(c) Neil Gershenfeld 1/23/26