Diffusion
Diffusion on a network, du_i/dt = D Σ_j A_ij (u_j - u_i), which spreads a unit mass to every node.
- Name:
network_diffusion - Field:
mathematics/networks - Formalism: ordinary differential equations
- References: Thurner et al. 2018; Jensen 2023
Factory
@register_model(
field='mathematics/networks',
tags=('network', 'diffusion', 'graph laplacian', 'ode'),
references=('Thurner et al. 2018', 'Jensen 2023'),
defaults={'until': 10.0, 'dt': 0.1},
)
def network_diffusion(graph: nx.Graph | None = None, *, D: float = 1.0, source: int = 0) -> System:
"""Return diffusion on a network, `du_i/dt = D Σ_j A_ij (u_j - u_i)`, which spreads a unit mass to every node.
The total is conserved, and on a connected graph the content of every node tends to `1 / n`.
Args:
graph:
The undirected NetworkX graph, or `None` for Zachary's karate club.
D:
The diffusion coefficient of every link.
source:
The index, in node order, of the node that holds the whole mass at the start.
Returns:
system:
The model, with the contents `nodes.u`.
Examples:
>>> import networkx as nx
>>> from skcomplex.models.mathematics.networks import network_diffusion
>>> network_diffusion(nx.path_graph(4)).elements[0].check().prepare().size_
4
"""
graph = _graph_or_default(graph)
initial = np.zeros(graph.number_of_nodes())
initial[source] = 1.0
nodes = FieldElement('nodes', element=NodeContentElement, space=NetworkSpace(graph), values={'u': initial, 'D': D})
diffusion = DiffusionInteraction('diffusion', target='nodes', state='u', coefficient=lambda s: s.target.D)
return System(elements=[nodes], interactions=[diffusion])
Components
class NodeContentElement(Element):
"""Content of a node of a network, which diffuses along the links with the coefficient `D`.
Args:
label:
The label of the node or the field.
visual:
How the element is drawn, or `None` for `default_visual_`.
u:
The initial content.
D:
The diffusion coefficient of every link.
Examples:
>>> from skcomplex.models.mathematics.networks import NodeContentElement
>>> NodeContentElement('node', D=0.5).D
0.5
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
u: Annotated[float, StateSpec()] = 0.0,
D: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
) -> None:
super().__init__(label, visual=visual)
self.u = u
self.D = D
class DiffusionInteraction(ConservationInteraction):
"""Diffusion `du/dt = D laplacian(u)` of a state `u` of a field, the conservation law of the flux `-D gradient(u)`.
The coefficient `D` is uniform in space, a number or a function of the namespace that returns an expression of
parameters and time. The interaction lowers to the Laplacian of the field, which its discretization evaluates, so
the total of the state on a periodic field stays constant and a point mass spreads as the heat kernel.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the field whose state diffuses.
state:
The name of the state of the field that diffuses.
coefficient:
A function of the namespace that returns the diffusion coefficient, or a number.
Raises:
ValidationError:
The state is not a name, or the coefficient is neither callable nor a number.
Examples:
>>> from typing import Annotated
>>> import numpy as np
>>> from skcomplex.core.elements import Element, FieldElement
>>> from skcomplex.core.interactions import DiffusionInteraction
>>> from skcomplex.core.spaces import LatticeSpace
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.specs import StateSpec
>>> class CellElement(Element):
... def __init__(self, label, u: Annotated[float, StateSpec()] = 0.0):
... super().__init__(label)
... self.u = u
>>> space = LatticeSpace((40,), spacing=0.1, boundary='periodic')
>>> field = FieldElement('rod', element=CellElement, space=space, values={'u': np.eye(40)[20]})
>>> heat = DiffusionInteraction('heat', target='rod', state='u', coefficient=0.1)
>>> trajectory = Simulation(System(elements=[field], interactions=[heat]), dt=0.01).check().prepare()
>>> round(float(trajectory.run(until=1.0).trajectory_.values('rod.u')[-1].sum()), 6)
1.0
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'u',
coefficient: Callable[[SimpleNamespace], object] | float = 1.0,
) -> None:
super().__init__(label, target=target, state=state)
self.coefficient = coefficient
def _check(self) -> None:
"""Raise an error for a state that is not a name, or a coefficient that is neither callable nor a number.
Raises:
ValidationError:
The type declares states, the state is not a name, or the coefficient is neither callable nor a number.
"""
super()._check()
if not (callable(self.coefficient) or isinstance(self.coefficient, int | float)):
error_msg = f'The coefficient of {self.label!r} is {self.coefficient!r}. Give a function or a number.'
raise ValidationError(error_msg)
def _get_transport_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the rate `D laplacian(u)` of the state, the divergence of the diffusive flux.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespace of the field.
Returns:
laws:
The rate law of the state.
"""
coefficient = self.coefficient(namespace) if callable(self.coefficient) else self.coefficient
state = getattr(self._get_member(namespace, 'target'), self.state)
return [RateLaw(state, coefficient * laplacian(state))]
Equations
nodes.u[i]: diff(u_nodes(t), t) == D_nodes*laplacian(u_nodes(t))