Differential equations

A rule for how something changes from moment to moment. scipy's solve_ivp follows it from a start you can drag.

Setup

Run this first.

import numpy as np
from functools import lru_cache
from scipy.integrate import solve_ivp

@lru_cache(maxsize=16)
def solve(f, start, t_end, steps=400):
    """The path from `start` over time 0 to t_end, as `steps` rows of state.
    Cached, so the path and the text above it solve once per drag, not twice.
    If the solver stops early, its last row repeats, so the path always has `steps` rows."""
    sol = solve_ivp(f, (0, t_end), start, t_eval=np.linspace(0, t_end, steps), rtol=1e-7, atol=1e-9)
    path = sol.y.T
    return np.vstack([path, np.repeat(path[-1:], steps - len(path), axis=0)])

def field(f, xs, ys, length=0.3, color='#78716c'):
    """At each grid point, a short line with a dot on the end that points the way f moves a point there."""
    X, Y = np.meshgrid(xs, ys)
    U, V = (np.broadcast_to(v, X.shape) for v in f(0, [X, Y]))
    size = np.hypot(U, V)
    size[size == 0] = 1
    tails = pie.point(X.ravel().tolist(), Y.ravel().tolist(), hidden=True)
    heads = pie.point((X + length * U / size).ravel().tolist(), (Y + length * V / size).ravel().tolist(), hidden=True)
    pie.line(tails, heads, stroke=color)
    pie.circle(heads, 0.05, fill=color, stroke='none')

def polyline(xs, ys, stroke):
    """Joins the points (xs[i], ys[i]) in order with straight lines."""
    pie.line(pie.point(xs[:-1].tolist(), ys[:-1].tolist(), hidden=True), pie.point(xs[1:].tolist(), ys[1:].tolist(), hidden=True), stroke=stroke)

A pendulum

x: angle. y: speed.

pie.clear()

def pendulum(t, s):
    angle, speed = s
    return [speed, -np.sin(angle)]

def pendulum_path(start):
    return solve(pendulum, start, 20)

def describe(start):
    angle, speed = start
    energy = speed ** 2 / 2 - np.cos(angle)
    if energy < 1:
        return f'swings up to {np.degrees(np.arccos(-energy)):.0f}° each side'
    return 'swings over the top and keeps spinning'

field(pendulum, np.arange(-4, 4.1, 0.5), np.arange(-3, 3.1, 0.5))
# The dividing line: a start on it just reaches the top, and stops there.
angles = np.linspace(-4, 4, 101)
for sign in (1, -1):
    polyline(angles, sign * 2 * np.abs(np.cos(angles / 2)), '#b45309')
start = pie.point(1, 0, fill='crimson')
path = pie.lift(pendulum_path, pie.Array[pie.Point])(start, hidden=True)
pie.line(path[:-1], path[1:], stroke='#15803d')
pie.text('angle', 4.4, -0.3, 12, stroke='#57534e')
pie.text('speed', 0, 3.4, 12, stroke='#57534e')
pie.text(pie.lift(describe, str)(start), 0, 4.8, 14, stroke='#1c1917')
pie.text('Drag the red point', 0, -4.8, 14, stroke='#1c1917')

Inside the amber line it swings. Outside, it spins over the top.

Predators and prey

x: rabbits. y: foxes. Foxes eat rabbits.

pie.clear()
pie.scalar(-3.5, out='originX')
pie.scalar(-3.5, out='originY')
a, b, c, d = 1, 0.4, 0.3, 0.9

def predator_prey(t, s):
    prey, predators = s
    return [a * prey - b * prey * predators, c * prey * predators - d * predators]

def predator_prey_path(start):
    return solve(predator_prey, start, 30)

def describe(start):
    if min(start) <= 0:
        return 'start with some of each'
    path = predator_prey_path(start)
    (prey_low, predators_low), (prey_high, predators_high) = path.min(axis=0), path.max(axis=0)
    return f'prey: {prey_low:.1f} to {prey_high:.1f}, predators: {predators_low:.1f} to {predators_high:.1f}'

field(predator_prey, np.arange(0.5, 8.1, 0.5), np.arange(0.5, 8.1, 0.5))
# The balance point: neither count changes.
pie.circle(pie.point(d / c, a / b, hidden=True), 0.15, stroke='#b45309')
start = pie.point(4.5, 2.5, fill='crimson')
path = pie.lift(predator_prey_path, pie.Array[pie.Point])(start, hidden=True)
pie.line(path[:-1], path[1:], stroke='#15803d')
pie.text('prey', 7.6, -0.3, 12, stroke='#57534e')
pie.text('predators', 0, 7.9, 12, stroke='#57534e')
pie.text(pie.lift(describe, str)(start), 3.5, 8.3, 14, stroke='#1c1917')
pie.text('Drag the red point', 3.5, -1.3, 14, stroke='#1c1917')

Every start loops around the amber ring, where neither count changes.

An orbit

A planet (red) around a sun (amber). The blue line is its starting speed.

pie.clear()
pie.scalar(0, out='gridOpacity')
GM = 4  # the strength of the sun's pull

def gravity(t, s):
    x, y, vx, vy = s
    r3 = (x * x + y * y + 0.01) ** 1.5  # + 0.01: no dividing by 0 at the sun's center
    return [vx, vy, -GM * x / r3, -GM * y / r3]

def state(planet, handle):
    """Position and speed: the speed is the blue line, from the planet to the handle."""
    return (planet[0], planet[1], handle[0] - planet[0], handle[1] - planet[1])

def orbit_path(planet, handle):
    return solve(gravity, state(planet, handle), 40, 800)[:, :2]

def describe(planet, handle):
    x, y, vx, vy = state(planet, handle)
    energy = (vx ** 2 + vy ** 2) / 2 - GM / np.hypot(x, y)
    spin = x * vy - y * vx
    stretch = np.sqrt(max(0, 1 + 2 * energy * spin ** 2 / GM ** 2))
    if spin ** 2 / (GM * (1 + stretch)) < 0.3:
        return 'falls into the sun'
    if energy >= 0:
        return 'escapes: too fast for the sun to hold'
    period = 2 * np.pi * np.sqrt((-GM / (2 * energy)) ** 3 / GM)
    return f'goes around every {period:.1f} time units' + (', nearly a circle' if stretch < 0.05 else '')

pie.circle(pie.point(0, 0, hidden=True), 0.3, fill='#b45309', stroke='none')
planet = pie.point(3, 0, fill='crimson')
handle = pie.point(3, 1, fill='#2563eb')
pie.line(planet, handle, stroke='#2563eb')
path = pie.lift(orbit_path, pie.Array[pie.Point])(planet, handle, hidden=True)
pie.line(path[:-1], path[1:], stroke='#15803d')
pie.text(pie.lift(describe, str)(planet, handle), 0, 4.8, 14, stroke='#1c1917')
pie.text('Drag the red or the blue point', 0, -4.8, 14, stroke='#1c1917')