Fractional Diffusion Solver

A pure-Python library for solving time-dependent fractional diffusion equations using eigenfunction expansion methods.

GitHub-Only Distribution | View on GitHub

⚠️ Warning

This software is experimental and under active development. It may contain bugs or undergo significant changes. Use with caution for research or production purposes.

Quick Start

Installation from GitHub

git clone https://github.com/example/fractional-diffusion.git
cd fractional-diffusion
pip install -e .

Optional Dependencies (for examples)

pip install -e ".[examples]"

What It Does

Solves the fractional diffusion equation: $$\frac{\partial u}{\partial t} + (-\Delta)^{\alpha} u = 0$$

with Dirichlet boundary conditions ($u(0,t) = u(1,t) = 0$) and initial condition $u(x, 0) = u_0(x)$.

Mathematical Formulation

The solution is obtained via eigenfunction expansion: $$u(x,t) = \sum_{n=1}^{\infty} a_n e^{-\mu_n t} \psi_n(x)$$

where:

• $\psi_n(x)$ are the eigenfunctions (e.g., $\sqrt{2}\sin(n\pi x)$ for the standard Laplacian)
• $\mu_n$ are the eigenvalues (roots of a characteristic polynomial depending on $\alpha$)
• $a_n = \int_0^1 u_0(x) \psi_n(x) , dx$ are the expansion coefficients

The key difference from simpler approaches is that $\mu_n$ are not simply $\lambda_n^{\alpha}$; instead, they satisfy a fractional eigenvalue problem. This distinction is crucial for accurate solutions, especially when $\mu_n$ is small.

Key Features

• Accurate eigenvalues: Uses $\mu_n$ (roots of a polynomial) rather than $\lambda_n^{\alpha}$—the difference is especially important for small $\mu_n$
• Flexible fractional powers: Supports any $\alpha > 0$
• Efficient computation: ~0.7ms per time point evaluation
• Boundary validation: Automatic checking of initial conditions

Example: Comparing Fractional Powers

Different fractional powers $\alpha$ produce dramatically different decay rates:

Core API

from fractional_diffusion import TimeEvolution
import numpy as np

# Setup: load or compute eigenvectors and eigenvalues
eigenvectors = ...  # shape (n_modes, n_points)
eigenvalues = ...   # shape (n_modes,)
x = np.linspace(0, 1, n_points)

# Create solver
solver = TimeEvolution(eigenvectors, eigenvalues, x)

# Solve with initial condition
initial_condition = np.sin(np.pi * x)
u_final = solver.solve(initial_condition, t=0.1)

# Or evaluate at multiple times
u_all, times = solver.solve(initial_condition, t=np.linspace(0, 1, 50), 
                             return_all_times=True)

Examples

The package includes 5 complete working examples:

1. Basic Usage (examples/01_basic_usage.py)

   • Single and multiple time evaluations
   • Simple API usage

2. Visualization (examples/02_visualization.py)

   • Step, Gaussian, triangular, multi-peak ICs
   • Real-time plotting of evolution

3. Advanced Features (examples/03_advanced.py)

   • Performance benchmarks (~0.7ms/eval)
   • Boundary condition checking
   • Memory efficiency analysis

4. Custom Eigendata (examples/04_custom_eigendata.py)

   • Integration with external eigensolvers
   • Fortran/FEM eigenvalue data
   • Orthonormalization guidance

5. Fractional Power Comparison (examples/05_fractional_powers_comparison.py)

   • Side-by-side comparison of α = 0.5, 1.0, 1.5
   • Temporal decay analysis
   • Quantitative results showing α-dependent behavior

Run Examples

cd fractional-diffusion
python examples/01_basic_usage.py
python examples/02_visualization.py
python examples/03_advanced.py
python examples/04_custom_eigendata.py
python examples/05_fractional_powers_comparison.py

Key Features

✓ Boundary Condition Validation

Automatic checking that initial conditions satisfy Dirichlet boundaries:

# Raises error if IC violates BC (configurable tolerance)
u = solver.solve(ic, t=0.1, check_bc=True, bc_tol=1e-6)

# Or disable checking for exploratory analysis
u = solver.solve(ic, t=0.1, check_bc=False)

✓ Batch Time Processing

Efficiently evaluate at multiple time points:

times = np.linspace(0, 1, 100)
u_all, t_eval = solver.solve(ic, t=times, return_all_times=True)
# u_all.shape = (100, n_points)

✓ Custom Eigensolvers

Works with any eigenvalue/eigenvector pair:

# From Fortran, FEM, or other solvers
evals, evecs = your_custom_eigensolver(L, bc)
solver = TimeEvolution(evecs, evals, x)

Mathematical Background

The solver uses eigenfunction expansion:

1. Expand IC: $a_n = \int_0^1 u_0(x) \psi_n(x) dx$
2. Solve ODE: $u_n(t) = a_n e^{-\mu_n^{\alpha} t}$
3. Reconstruct: $u(x,t) = \sum_n u_n(t) \psi_n(x)$

Fractional powers scale eigenvalues: $$\mu_n^{\alpha} \approx \lambda_n^{(\alpha)} = (\lambda_n)^{\alpha}$$ with $\lambda_n = (n \pi)^2$ Result: Decay rate directly depends on fractional power $\alpha$

Installation Requirements

• Python ≥ 3.8
• NumPy ≥ 1.20
• SciPy ≥ 1.6
• Matplotlib ≥ 3.5 (optional, for visualization examples)

Project Structure

fractional-diffusion/
├── fractional_diffusion/       # Main package
│   ├── __init__.py
│   └── time_evolution.py       # TimeEvolution class
├── examples/                   # 5 working examples
│   ├── 01_basic_usage.py
│   ├── 02_visualization.py
│   ├── 03_advanced.py
│   ├── 04_custom_eigendata.py
│   └── 05_fractional_powers_comparison.py
├── pyproject.toml              # Package metadata
├── README_GITHUB.md            # This file
└── LICENSE                     # MIT License

References

• Kundu, A., Bernardin, C., Saito, K., Kundu, A., & Dhar, A. (2019). Fractional equation description of an open anomalous heat conduction set-up. Journal of Statistical Mechanics: Theory and Experiment, 2019(1), 013205. https://doi.org/10.1088/1742-5468/aaf630
• Kwaśnicki, M. Ten Equivalent Definitions of the Fractional Laplace Operator. FCAA 20, 7–51 (2017). https://doi.org/10.1515/fca-2017-0002

Citation

If you use this software in research, please cite:

@software{fractional_diffusion_github,
  title={Fractional Diffusion: Time Evolution Solver},
  author={Aritra Kundu},
  year={2017-2026},
  howpublished={\url{https://github.com/akunduphys/FractionalLaplacian}}
}

License

MIT License - see LICENSE for details

Support

For issues, questions, or contributions, visit the GitHub repository

---

Built with Python 3.8+ • Pure NumPy/SciPy • MIT Licensed

Luxembourg, 2026