Mathematical Computing

The Nature of Mathematical Modeling (draft)

Open- vs closed-source

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: pip, conda, uv, ...
    • pip install jupyterlab, jupyterlab-git, jupyterlite-core, jupyterlite-pyodide-kernel, ipympl, ipywidgets, ipycanvas, nbconvert, pretty-jupyter, numpy, scipy, sympy, numba, pandas, jax, scikit-learn, matplotlib, plotly, pythreejs, anywidget, opencv-contrib-python

Rust

  • Fixes recurring language security and stability issues: buffer overflows, memory leaks, race conditions, ...
  • Compiled, can target multiple platforms
  • tutorial

Web

Javascript/ECMAScript

WebAssembly (WASM)

WebGPU

  • Graphics and compute acceleration

Coding

Text editors

  • small, fast, can eliminate need to use mouse
  • Vi
  • Emacs
  • ...

Integrated Development Environments (IDEs)

  • adds support for documentation, development, debugging
  • VS Code
  • Eclipse
  • ...

Large Language Models

  • enables vibe coding, describing rather than writing code; can make mistakes, infringe copyrights, ...
  • Cursor
  • Codex
  • Copilot
  • ...

Notebooks

  • Combine coding, development, documentation, distribution, collaboration

Jupyter

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

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

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

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.591733, estimated MFlops = 84.497606

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.065493, estimated MFlops = 7634.419014

Numba parallel version
NPTS = 10000000000, pi = 3.141593
time = 0.570371, estimated MFlops = 87662.181046

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.060820, estimated MFlops = 822.103051

Jax

  • Python just-in-time compilation, parallel acceleration, automatic differentiation
print('Jax pi calculation ...\n')
#
import os
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.127273
run Jax version:
NPTS = 100000000, pi = 3.141592
time = 0.056738, estimated MFlops = 8812.451676

compile Jax parallel version:
time = 0.366621
run Jax parallel version using 16 CPUs:
NPTS = 1600000000, pi = 3.141592
time = 0.411482, estimated MFlops = 19441.937702

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
import numpy as np
import sympy as sy
print('sqrt(8):')
print(f"NumPy: {np.sqrt(8)}")
print(f"SymPy: {sy.sqrt(8)}")

sqrt(8):
NumPy: 2.8284271247461903
SymPy: 2*sqrt(2)

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
x = np.arange(-l,l,0.2)
y = np.sin(k*x)/(k*x)
plt.ion()
fig,ax = plt.subplots()
fig.canvas.header_visible = False
line, = ax.plot(x,y)
plt.title('sin(kx)/(kx)')
plt.show()
#
# handle slider changes
#
def slider_handler(change):
    k = change['new']
    y = np.sin(k*x)/(k*x)
    line.set_ydata(y)
    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)
Figure 
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': '6eb63721-16c2-4b47-a507-34851cc3fc67',
              '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…

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

  1. Simulate/plot/render/... bouncing balls

(c) Neil Gershenfeld 1/4/26