asymptotics

A symbolic perturbation theory toolkit for Python

Write your perturbation problem as a string. Get symbolic order-by-order results,
LaTeX display, numerical verification, and direct evaluation — automatically.

What is asymptotics?

asymptotics is a Python library that automates classical perturbation methods for algebraic equations and ordinary differential equations (ODEs). Instead of manually deriving order-by-order equations, substituting ansätze, and solving each level — you write your problem as a plain string and call a method:

from asymptotics import ODE

eq  = ODE("u'' + u + eps*u**3", small_param="eps",
          conditions=["u(0) = 1", "u'(0) = 0"])

sol = eq.expand_lindstedt(order=2)
sol.show()       # beautiful LaTeX in Jupyter; clean text in terminal

The library handles the rest: symbolic expansion, secular term elimination, condition application, and optional numerical comparison against SciPy.

Unlike black-box tools (Mathematica's AsymptoticDSolveValue), asymptotics exposes every intermediate step — each order's ODE, its general and particular solution, secular terms, frequency corrections — as live SymPy expressions you can inspect and manipulate.

Installation

pip install asymptotics

For local development with example notebooks:

git clone https://github.com/saadgroup/asymptotics
cd asymptotics
pip install -e ".[dev,notebook]"

Requirements: Python ≥ 3.10 · SymPy ≥ 1.12 · NumPy ≥ 1.24 · SciPy · Matplotlib

Methods at a glance

Class	Method	Use when…
`AlgebraicEquation`	`.expand_regular(order)`	Nonlinear algebraic equation $f(x,\varepsilon)=0$
`AlgebraicSystem`	`.expand_regular(order)`	Coupled algebraic system
`ODE`	`.expand_regular(order)`	ODE with small nonlinear/forcing term (IVP or BVP)
`ODE`	`.expand_lindstedt(order)`	Nonlinear oscillator — strains time to remove secular terms
`ODE`	`.expand_multiple_scales(order)`	Oscillator with slow amplitude/phase modulation
`ODE`	`.expand_boundary_layer()`	Singular BVP — $\varepsilon$ multiplies highest derivative
`ODE`	`.begin_expansion(order)`	Step-by-step control when SymPy cannot solve a leading order
`ODESystem`	`.expand_regular(order)`	Coupled system of ODEs

Every result supports a consistent four-method API:

Method	What it does
`sol.show()`	LaTeX display in Jupyter; plain text in terminal
`sol.eval(eps, at=...)`	Evaluate composite as a NumPy array
`sol.compare_numeric(eps)`	Numerical verification plot via SciPy
`sol.to_latex(...)`	Export ready-to-paste LaTeX source

Quick examples

Algebraic equation

from asymptotics import AlgebraicEquation

eq  = AlgebraicEquation("x**3 + eps*x - 1", dependent="x", small_param="eps")
sol = eq.expand_regular(order=3)
sol.show()
# x(ε) = 1 - ε/3 - ε³/81 + O(ε⁴)

sol[0].solution    # x₀ = 1
sol[1].solution    # x₁ = -1/3
sol.composite      # full SymPy expression

sol.eval(eps=0.1)                          # float
sol.eval(eps=[0.1, 0.2, 0.3])             # ndarray
sol.compare_numeric(eps=0.3)              # plot vs scipy.fsolve

Coupled algebraic system

from asymptotics import AlgebraicSystem

sys = AlgebraicSystem(
    equations   = ["x**2 + eps*sin(y) - 1", "y - eps*cos(x)"],
    dependents  = ["x", "y"],
    small_param = "eps",
)
sol = sys.expand_regular(order=3)
sol.show()

sol["x"].composite   # expansion for x(ε)
sol["y"].composite   # expansion for y(ε)

Regular perturbation — IVP

from asymptotics import ODE

# dependent='u' and independent='t' inferred automatically from conditions
eq = ODE(
    "u'' + u + eps*u**3",
    small_param = "eps",
    conditions  = ["u(0) = 1", "u'(0) = 0"],
)
sol = eq.expand_regular(order=2)
sol.show()

sol[1].secular   # True — secular terms detected; use Lindstedt or multiple scales

Regular perturbation — BVP

# BVP detected automatically from two distinct boundary points
# independent='x' inferred
eq = ODE(
    "u'' + eps*u",
    small_param = "eps",
    conditions  = ["u(0) = 0", "u(1) = 1"],
)
sol = eq.expand_regular(order=3)
sol.show()

Lindstedt–Poincaré

Fix secular terms in nonlinear oscillators by straining the time coordinate $\tau = \omega(\varepsilon),t$.

eq = ODE(
    "u'' + u + eps*u**3",      # Duffing oscillator
    small_param = "eps",
    conditions  = ["u(0) = 1", "u'(0) = 0"],
)
sol = eq.expand_lindstedt(order=2)
sol.show()

sol.omega_0           # ω₀ = 1   (auto-detected from unperturbed equation)
sol.omega_expansion   # ω(ε) = 1 + 3ε/8 - 21ε²/256
sol[1].omega_k_val    # ω₁ = 3/8  (classical result)
sol.composite_t       # u(t, ε) — uniformly valid in physical time t

import numpy as np
t_vals = np.linspace(0, 40, 500)
sol.eval(eps=0.1, at=t_vals)        # ndarray
sol.compare_numeric(eps=0.1)        # plot vs scipy.solve_ivp

Works for any natural frequency — detected automatically from the unperturbed equation:

eq = ODE("u'' + 4*u + eps*u**3", small_param="eps",
         conditions=["u(0) = 1", "u'(0) = 0"])
sol = eq.expand_lindstedt(order=2)
sol.omega_0   # 2 — detected from u'' + 4u = 0

Multiple scales

For problems where amplitude or phase evolve slowly.

# Weakly damped oscillator
eq = ODE(
    "u'' + u + eps*u'",
    small_param = "eps",
    conditions  = ["u(0) = 1", "u'(0) = 0"],
)
sol = eq.expand_multiple_scales(order=1)
sol.show()

sol.amplitude_A   # A(T₁) = e^{-T₁/2}  — exact Bernoulli solution
sol.composite_t   # e^{-εt/2} · cos(t)  — matches exact result at O(1)

# Van der Pol oscillator — limit cycle
eq = ODE(
    "u'' + u + eps*(u**2 - 1)*u'",
    small_param = "eps",
    conditions  = ["u(0) = 1", "u'(0) = 0"],
)
sol = eq.expand_multiple_scales(order=1)
sol.amplitude_A   # 2√(eᵀ¹/(eᵀ¹+3))  →  2  as T₁→∞  (limit cycle)

Boundary layers

For singular BVPs where $\varepsilon$ multiplies the highest derivative. The layer location is detected automatically from the sign of $p(x)$.

eq = ODE(
    "eps*u'' + u' + u",   # p(0) = 1 > 0  →  layer at x = 0
    small_param = "eps",
    conditions  = ["u(0) = 0", "u(1) = 1"],
)
sol = eq.expand_boundary_layer()
sol.show()

sol.layer_location   # 'x = 0'
sol.outer            # outer solution (away from layer)
sol.inner            # inner solution U(ξ) in stretched coord ξ = x/ε
sol.composite        # u_out + u_in − u_match  (Van Dyke rule)

sol.compare_numeric(eps=0.05)

Variable coefficients are fully supported:

eq = ODE(
    "eps*u'' + (1 + x)*u' - u",
    small_param = "eps",
    conditions  = ["u(0) = 1", "u(1) = 2"],
)

Coupled ODE system

from asymptotics import ODESystem

sys = ODESystem(
    equations   = ["u' + u + eps*v", "v' + 2*v + eps*u**2"],
    dependents  = ["u", "v"],
    small_param = "eps",
    independent = "t",
    conditions  = ["u(0) = 1", "v(0) = 1"],
)
sol = sys.expand_regular(order=2)
sol.show()

sol["u"].composite              # full expansion for u(t, ε)
sol["v"].composite              # full expansion for v(t, ε)
sol["u"][1].particular_solution # u₁(t)

import numpy as np
t_vals = np.linspace(0, 10, 200)
sol.eval(eps=0.1, at=t_vals)    # {'u': ndarray, 'v': ndarray}
sol.compare_numeric(eps=0.1)

Works for any number of coupled equations.

Step-by-step API

When SymPy cannot solve a leading-order equation (e.g. a nonlinear ODE at O(1)), begin_expansion() gives you full manual control while letting the library handle all linear higher-order equations automatically.

sol = eq.begin_expansion(order=2)

# Inspect all equations immediately — nothing is solved yet
sol.show()

# Examine order-k ODE
sol[0].ode    # O(1) equation — purely symbolic
sol[1].ode    # O(ε) equation — symbolic if O(1) unsolved;
              #                  fully substituted if O(1) is solved

# Solve or provide solutions
sol[0].solve()                      # try SymPy — fails gracefully with clear message
sol[0].set_solution("4*eta*(1-eta)") # provide expression as string or SymPy expr
sol[1].solve()                      # SymPy handles linear higher orders
sol.solve_all()                     # attempt all remaining

# Inspect state
sol[k].is_solved                    # bool
sol[k].particular_solution          # SymPy expression (if solved)
sol[k].general_solution             # SymPy expression (if solved)
sol.n_solved                        # number of solved orders
sol.n_pending                       # number of unsolved orders

# Full standard API once all orders are solved
sol.composite
sol.show()
sol.to_latex()
sol.eval(eps=0.1, at=t_vals)
sol.compare_numeric(eps=0.1)

If composite, eval, compare_numeric, or to_latex is called before all orders are solved, a clear NotReadyError lists exactly which orders are pending.

Non-standard gauge sequences

By default the ansatz uses the standard power sequence ${1, \varepsilon, \varepsilon^2, \ldots}$. You can supply a non-standard gauge when the problem calls for it:

# Half-power gauge: {1, √ε, ε, ε^(3/2), ...}
sol = eq.expand_regular(order=3, gauge="sqrt(eps)")

# Logarithmic gauge: {1, log(ε), log²(ε), ...}
sol = eq.expand_regular(order=2, gauge=["1", "log(eps)", "log(eps)**2"])

# Inspect the gauge used
sol.gauge   # list of SymPy gauge functions

Symbolic parameters

Parameters other than the dependent variable, independent variable, and small parameter are detected automatically and a warning is issued:

eq = ODE("u'' + a*u + eps*u**3", small_param="eps",
         conditions=["u(0) = 1", "u'(0) = 0"])
# ⚠️  symbolic parameters detected: {'a'}

# Parameters also work in conditions
eq = ODE("u'' + u + eps*u**3", small_param="eps",
         conditions=["u(0) = A", "u'(0) = 0"])
# ⚠️  symbolic parameters detected: {'A'}

# show() and to_latex() always work — results stay symbolic
sol.show()

# Provide values at eval/compare time
sol.eval(eps=0.1, at=t_vals, params={"a": 2.0, "A": 1.5})
sol.compare_numeric(eps=0.1, params={"a": 2.0, "A": 1.5})
# Missing params → clear error listing all required names

The four core methods

`sol.show()`

LaTeX rendering in Jupyter; clean plain text in the terminal. The small parameter is always displayed as $\varepsilon$, regardless of the name you chose (eps, mu, delta, …).

sol.show()                   # full hierarchy
sol.show(orders=[0, 1])      # selected orders only
sol.show(mode='text')        # force plain text (useful in scripts)

`sol.eval(eps, at=None, params=None)`

Evaluate the composite expansion and return a NumPy array.

import numpy as np
t_vals = np.linspace(0, 20, 400)

# ODE — single eps
u = sol.eval(eps=0.1, at=t_vals)             # ndarray

# ODE — multiple eps values
u = sol.eval(eps=[0.1, 0.2], at=t_vals)      # dict {0.1: ndarray, 0.2: ndarray}

# Algebraic
x = sol.eval(eps=0.1)                        # float
x = sol.eval(eps=[0.1, 0.2, 0.3])           # ndarray

# ODESystem
r = sol.eval(eps=0.1, at=t_vals)            # {'u': ndarray, 'v': ndarray}

# With symbolic parameters
sol.eval(eps=0.1, at=t_vals, params={"a": 2.0})

For Lindstedt and multiple scales, composite_t (solution in physical time $t$) is used automatically — no need to handle the strained-time symbol yourself.

`sol.compare_numeric(eps, params=None, plot_range=None, filename=None)`

Generate a comparison plot between the perturbation expansion and a SciPy numerical reference solution.

sol.compare_numeric(eps=0.1)
sol.compare_numeric(eps=[0.1, 0.2, 0.3])           # overlay multiple eps
sol.compare_numeric(eps=0.1, params={"a": 2.0})
sol.compare_numeric(eps=0.1, plot_range=[0, 40])    # override default range
sol.compare_numeric(eps=0.1, filename="fig.pdf")    # save to file

The plot range is inferred automatically from the problem's conditions. No problem= argument is required — the problem is stored on the hierarchy.

Returns a dict with keys that depend on problem type:

Problem type	Keys
All ODE methods	`'t'`, `'u_pert'`, `'u_numerical'`, `'fig'`
Boundary layer	additionally `'u_outer'`, `'u_inner'`, `'u_composite'`
ODE system	`'u_pert'` and `'u_numerical'` are dicts keyed by variable name
Algebraic	`'eps'`, `'perturbation'`, `'numerical'`, `'fig'`

SciPy solvers used: solve_ivp (IVP, Lindstedt, multiple scales), solve_bvp (BVP, boundary layer), fsolve (algebraic).

`sol.to_latex(environment='align', show_orders=False, filename=None)`

Export results as ready-to-paste LaTeX source. The small parameter is always rendered as \varepsilon.

sol.to_latex()                              # print to console
sol.to_latex(filename="result.tex")        # save to file
sol.to_latex(show_orders=True)             # include u₀, u₁, u₂, … before composite
sol.to_latex(environment='equation')       # default: 'align'
sol.to_latex(environment='gather')

Lindstedt output includes the frequency expansion $\omega(\varepsilon)$, the composite in strained time $\tau$, and the composite in physical time $t$. Boundary layer output includes outer, inner, and composite separately.

Inspecting intermediate steps

Every hierarchy exposes the full symbolic work at each order:

sol = eq.expand_regular(order=3)

# Per-order entries
sol[k].ode                   # the ODE/equation at order k
sol[k].general_solution      # with free integration constants
sol[k].particular_solution   # constants fixed by ICs/BCs
sol[k].secular               # True if secular terms detected (ODE only)

# Assembly
sol.composite                # assembled expansion u(t, ε) as SymPy expression
sol.small_param              # the ε symbol

Lindstedt-specific attributes:

sol.omega_0                  # unperturbed frequency (auto-detected)
sol.omega_expansion          # ω(ε) as a SymPy series
sol[k].omega_k_val           # frequency correction ωₖ at order k
sol[k].secularity_condition  # the equation that determined ωₖ
sol.composite_t              # u(t, ε) with τ = ω(ε)·t substituted

Multiple scales-specific attributes:

sol.T0, sol.T1               # fast and slow time symbols
sol.amplitude_A              # A(T₁) — solved if B = 0 and Bernoulli ODE applies
sol.amplitude_B              # B(T₁)
sol[k].solvability_A         # amplitude ODE for A

Condition syntax

Conditions are plain strings — the same notation you'd write on paper:

conditions = ["u(0) = 1"]                          # 1st-order IVP
conditions = ["u(0) = 1", "u'(0) = 0"]             # 2nd-order IVP
conditions = ["u(0) = 0", "u(1) = 1"]              # BVP
conditions = ["u(pi) = 0", "u'(0) = sqrt(2)"]      # symbolic points/values
conditions = ["u(0) = 1/2", "u'(0) = 0"]           # rational values
conditions = ["u(0) = 0.9", "u'(0) = 0"]           # floats auto-converted

asymptotics automatically:

Detects IVP vs BVP from the number of distinct boundary points
Infers the dependent variable name from condition strings
Infers the independent variable (t for IVPs, x for BVPs) from the equation
Reports exactly what was inferred (with override syntax) at construction time
Raises clear ConditionError for wrong count, conflicts, or inconsistencies

Limit conditions for regularity at singular points:

conditions = [
    "F(0) = 0",
    "F(1) = 1",
    "lim(sqrt(2*eta)*F'', eta, 0) = 0",   # lim(expr, var, point) = value
]

At each order the library substitutes the general solution, identifies which free constants cause divergence, and sets them to zero automatically.

Auto-inference of variables

For ODE, dependent and independent are both optional:

# Fully minimal — both inferred
eq = ODE("u'' + u + eps*u**3", small_param="eps",
         conditions=["u(0) = 1", "u'(0) = 0"])
# ℹ️  dependent = 'u' (inferred from conditions)
# ℹ️  independent = 't' (IVP default)
# To override: ODE(..., dependent='u', independent='t')

# Override when the equation contains a non-standard independent variable
eq = ODE("u'' + sin(tau)*u + eps*u**3", small_param="eps",
         conditions=["u(0) = 1", "u'(0) = 0"],
         independent="tau")

Independent variable candidates inferred from the equation: {x, y, z, t}. All other symbols are treated as parameters. If a candidate symbol is ambiguous (could be independent variable or parameter), a hard error asks you to specify independent explicitly.

Higher-order ODEs

ODEs up to 6th order are supported using prime notation:

# 4th-order ODE
eq = ODE(
    "eta*F'''' + alpha*(eta*F''' + 2*F'') + eps/2*(F*F''' - F'*F'') + 2*F'''",
    small_param = "eps",
    independent = "eta",
    conditions  = ["F(0) = 0", "F(1/2) = 1", "F'(1/2) = 0",
                   "lim(sqrt(2*eta)*F'', eta, 0) = 0"],
)
sol = eq.begin_expansion(order=1)

Prime notation: u' (1st), u'' (2nd), u''' (3rd), u'''' (4th), u''''' (5th), u'''''' (6th).

Lindstedt and multiple scales are limited to 2nd-order oscillators. Regular perturbation and begin_expansion work for any supported order.

Error messages

asymptotics raises clear, actionable errors:

ConditionError: 3 conditions provided but the ODE is order 2 — need exactly 2.
  IVP: ["u(0) = 1", "u'(0) = 0"]
  BVP: ["u(0) = 0", "u(1) = 1"]

ConditionError: Conflicting conditions at t=0:
  u(0) = 1
  u(0) = 2

ValueError: Could not parse equation: 'u^2 + eps*u - 1'
  Use ** for powers:  u**2  not  u^2

NotReadyError: Cannot access 'composite' — 1 order(s) not yet solved: [0]
  Call sol[0].solve() or sol[0].set_solution(expr) first.

Project layout

asymptotics/
├── __init__.py                      ← public API
├── numerics.py                      ← compare_numeric() dispatcher + SciPy solvers
├── eval.py                          ← eval() for all hierarchy types
├── latex_export.py                  ← to_latex() for all hierarchy types
├── gauge.py                         ← non-standard gauge sequence support
├── core/
│   ├── problem.py                   ← AlgebraicEquation, AlgebraicSystem, ODE
│   ├── ode_system.py                ← ODESystem
│   ├── hierarchy.py                 ← OrderHierarchy, OrderEntry
│   ├── system_hierarchy.py          ← SystemHierarchy (coupled algebraic)
│   ├── conditions.py                ← parser: ParsedCondition, LimitCondition
│   └── exceptions.py               ← custom exceptions
├── methods/
│   ├── regular_algebraic.py         ← AlgebraicEquation solver
│   ├── regular_algebraic_system.py  ← AlgebraicSystem solver
│   ├── regular_ode.py               ← ODE regular perturbation
│   ├── regular_ode_system.py        ← ODESystem solver
│   ├── lindstedt.py                 ← Lindstedt–Poincaré
│   ├── multiple_scales.py           ← Method of multiple scales
│   ├── boundary_layer.py            ← Matched asymptotic expansions
│   └── stepwise.py                  ← StepwiseHierarchy, begin_expansion
├── display/
│   ├── jupyter.py                   ← algebraic LaTeX display
│   ├── ode_display.py
│   ├── ode_system_display.py
│   ├── lindstedt_display.py
│   ├── multiple_scales_display.py
│   └── boundary_layer_display.py
└── tests/                           ← 195 tests, all passing
    ├── test_regular_algebraic.py
    ├── test_algebraic_system.py
    ├── test_errors.py
    ├── test_gauge.py
    ├── test_ode.py
    ├── test_lindstedt.py
    ├── test_multiple_scales.py
    ├── test_boundary_layer.py
    └── test_ode_system.py

Example notebooks

Ten Jupyter notebooks covering every feature with worked examples:

Notebook	Topic
`01_introduction.ipynb`	Algebraic equations — root selection, symbolic parameters
`02_transcendental.ipynb`	Transcendental equations and coupled algebraic systems
`03_ode.ipynb`	Regular perturbation for ODEs — IVP and BVP, secular detection
`04_lindstedt.ipynb`	Lindstedt–Poincaré — Duffing oscillator, amplitude dependence
`05_multiple_scales.ipynb`	Multiple scales — damped oscillator, Van der Pol limit cycle
`06_boundary_layers.ipynb`	Matched asymptotic expansions — left/right layers, variable coefficients
`07_ode_system.ipynb`	Coupled ODE systems, predator-prey
`08_advanced.ipynb`	Advanced features and edge cases
`08_stepwise.ipynb`	Step-by-step API, higher-order ODEs, limit BCs
`09_gauge.ipynb`	Non-standard gauge sequences

Running the tests

pytest                      # all 195 tests
pytest -v                   # verbose output
pytest --cov=asymptotics    # with coverage report

Design philosophy

String-based input — write "u'' + u + eps*u**3", not u.diff(t, 2) + u + eps*u**3
Transparent by default — every intermediate step is a live SymPy expression you can inspect and reuse
Consistent API — expand_*, show(), eval(), compare_numeric(), to_latex() work identically across all problem types
Self-contained results — the problem is stored on the hierarchy; no need to pass it again to compare_numeric or eval
Fail clearly — errors name the problem, quote your input, and suggest a fix
No SymPy wrestling — the library manages all symbolic machinery internally; users stay at the mathematical level

Citation

If you use asymptotics in published work, please cite:

Tony Saad, asymptotics: A symbolic perturbation theory toolkit for Python, University of Utah, 2025. https://github.com/saadgroup/asymptotics

Name		Name	Last commit message	Last commit date
Latest commit History 19 Commits
asymptotics		asymptotics
examples		examples
.gitignore		.gitignore
LICENSE		LICENSE
README.md		README.md
pyproject.toml		pyproject.toml

Folders and files

Latest commit

History

Repository files navigation

asymptotics

What is asymptotics?

Installation

Methods at a glance

Quick examples

Algebraic equation

Coupled algebraic system

Regular perturbation — IVP

Regular perturbation — BVP

Lindstedt–Poincaré

Multiple scales

Boundary layers

Coupled ODE system

Step-by-step API

Non-standard gauge sequences

Symbolic parameters

The four core methods

sol.show()

sol.eval(eps, at=None, params=None)

sol.compare_numeric(eps, params=None, plot_range=None, filename=None)

sol.to_latex(environment='align', show_orders=False, filename=None)

Inspecting intermediate steps

Condition syntax

Auto-inference of variables

Higher-order ODEs

Error messages

Project layout

Example notebooks

Running the tests

Design philosophy

Citation

License

About

Resources

License

Uh oh!

Stars

Watchers

Forks

Releases

Packages 0

Uh oh!

Contributors

Uh oh!

Languages

`sol.show()`

`sol.eval(eps, at=None, params=None)`

`sol.compare_numeric(eps, params=None, plot_range=None, filename=None)`

`sol.to_latex(environment='align', show_orders=False, filename=None)`

Packages