Demathpy

A Python library for parsing and safely evaluating symbolic Ordinary and Partial Differential Equations (ODEs/PDEs) on numerical grids.

This repository provides:

• A class-based PDE solver (from demathpy import PDE) for running simulations.
• A lightweight symbol normalizer that converts human-readable mathematical notation into valid Python expressions.
• Equation-based Boundary Conditions and Initial Conditions.
• Built-in support for common differential operators and vector calculus notation using Finite Differences.

Key Features

1. The PDE Class

The core of the library is the PDE class. It manages the grid state, parsing, and time-stepping.

from demathpy import PDE
import numpy as np

# Define a Heat Equation: du/dt = Laplacian(u)
p = PDE("du/dt = lap(u)", space_axis=["x", "y"])

# Configure the grid
p.init_grid(width=100, height=100, dx=0.5)

# Set Initial Conditions using equations
p.initial = ["u = exp(-(x-25)**2 - (y-25)**2)"]
p.set_initial_state()

# Set Boundary Conditions using equations
# Format: "axis(coord) = value" or "axis=coord = value"
p.boundry = [
    "x=0 = 1.0",     # Left boundary (x=0) is fixed at 1.0
    "x=100 = 0.0",   # Right boundary (x=100) is fixed at 0.0
    "y=0 = 0.0",     # Bottom boundary is 0.0
    "periodic"       # Other unset boundaries (y=100) default to periodic or 0
]

# Run simulation
for _ in range(100):
    p.step(dt=0.01)

print(p.u.mean())

2. Equation-Based Configuration

You can configure boundaries and initial states using string equations instead of manual array manipulation.

Boundary Conditions (p.boundry list):

• Dirichlet: x=0 = 1.0 (Fixes value at boundary)
• Neumann: dx(u) = 0 (Not fully exposed yet, currently defaults to Dirichlet logic if value provided).
• Periodic: Use periodic keyword or leave empty for default periodic behavior (if implemented).

Initial Conditions (p.initial list):

• Scalar: u = sin(x) * cos(y)
• Vector Components: ux = 1.0, uy = 0.0 (if p.u_shape = ["ux", "uy"])

3. Vector Fields

The solver supports vector-valued PDEs (e.g., Navier-Stokes, Reaction-Diffusion systems).

# 2D Advection: du/dt = - (u · ∇) u
p = PDE("du/dt = -advect(u, u)", space_axis=["x", "y"])
p.u_shape = ["ux", "uy"] # Define component names
p.init_grid(width=50, height=50, dx=1.0)

# Initialize Vortex
p.initial = [
    "ux = -sin(y)",
    "uy = sin(x)"
]
p.set_initial_state()

Supported Operators

The parser recognizes and maps these to NumPy finite difference functions:

• Derivatives: du/dt, dx(u), dy(u), dz(u)
• Second Derivatives: dxx(u), dzz(u)
• Laplacian: lap(u) or ∇²u
• Gradient: grad(u) or ∇u (Returns vector)
• Divergence: div(u) or ∇·u (Expects vector input)
• Advection: advect(velocity, field) -> (velocity · ∇) field
• Math Functions: sin, cos, exp, log, abs, sqrt, tanh ...

Symbol Normalization

The parser supports Unicode and mathematical shorthand:

• α, β, γ → alpha, beta, gamma
• u² → u**2
• |u| → abs(u)

Workflow & Visualization

To integrate Demathpy into visualization software or interactive notebooks, you can use the get_grid() method to probe the field dynamics without advancing the simulation time.

Visualization Step-by-Step

1. Initialize:

   p = PDE("du/dt = lap(u) - u**3 + u", space_axis=["x", "y"])
   p.init_grid(width=20, height=20, dx=0.5)
   p.initial = ["u = 0.1 * sin(x)"] 
   p.boundry = ["periodic"]
   p.set_initial_state()
   

2. Probe the Vector Field (du/dt): Use get_grid(dt=0) to get the instantaneous rate of change. This is useful for visualizing flow fields or phase plots.

   # Get Rate of Change (RHS of PDE)
   du_dt = p.get_grid(dt=0)
   
   # Or calculate the hypothetical next step delta
   delta_u = p.get_grid(dt=0.01)
   

3. Predict on Arbitrary States: You can evaluate the PDE on a hypothetical state u_test without updating the solver's internal state. This is useful for drawing vector fields in phase space.

   # Create a test state
   test_u = np.sin(p.u) 
   
   # Calculate how the PDE would evolve this state
   # Returns the rate of change for the test state
   response = p.get_grid(u_state=test_u, dt=0)
   

4. Run Simulation Loops:

   import matplotlib.pyplot as plt
   
   for i in range(100):
       p.step(dt=0.01)
       if i % 10 == 0:
           plt.imshow(p.u) # Visualization logic
           # plt.show()
   

The ODE Class

Demathpy also efficiently solves Ordinary Differential Equations (ODEs) where the state depends only on time $t$. The API is identical to the PDE class.

from demathpy import ODE

# 1. EXPONENTIAL DECAY
# dy/dt = -y
o = ODE("dy/dt = -y", u_shape=["y"])
o.initial = ["y = 1.0"] 
o.init_state(shape=(1,)) # scalar system
o.set_initial_state()

for _ in range(100):
    o.step(dt=0.01) 
print(o.u) # Should be close to exp(-1)

# 2. VECTOR SYSTEMS (Predator-Prey)
# du/dt = u - u*v
# dv/dt = u*v - v
pp = ODE("du/dt = [u[0] - u[0]*u[1], u[0]*u[1] - u[1]]", u_shape=["u", "v"])
pp.initial = ["u = 1.1", "v = 1.0"]
pp.init_state(shape=(1,)) # Single ecosystem
pp.set_initial_state()

pp.step(dt=0.1)

# 3. BATCHED EXECUTION
# Simulating 1000 identical particles with different initial conditions
particles = ODE("dx/dt = -x + noise") # noise not impl by default but external vars work
# Or just decay
batch = ODE("dy/dt = -y")
batch.init_state(shape=(1000,)) # 1000 systems
# Set random initial states directly (or use equation if supported)
import numpy as np
batch.u[:] = np.random.rand(1000)

batch.step(0.1)

Key ODE Features:

• Equation Parsing: Supports dy/dt, d^2y/dt^2, vector syntax [a, b].
• Initialization: Use init_state(shape=...) where shape defines the batch size (independent systems).
• Probing: get_grid(u_state=..., dt=0) works exactly like PDE for generating phase portraits (return vector field at state).
• Functions: Includes sin, cos, exp, step, heaviside, sign, abs ...

License

MIT