Experimental formal proof
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
def dynamic_entropy_model(t, p, W, alpha):
“”"
ODE function for the 3-state DEM with a simple feedback control.
Parameters
----------
t : float
Current time (not used here explicitly, but required by solve_ivp).
p : array_like of shape (3,)
Current probability distribution [p1, p2, p3].
W : 2D array of shape (3, 3)
Transition rate matrix, where W[i,j] is rate from i -> j.
Diagonal entries can be zero or determined by row sums.
alpha : float
Feedback gain.
Returns
-------
dpdt : np.ndarray of shape (3,)
Time derivative of p = [dp1/dt, dp2/dt, dp3/dt].
"""
# Convert p to a numpy array for safety
p = np.array(p)
# Master-equation term (without control)
# dp_i/dt = sum_j [W[j,i]*p_j - W[i,j]*p_i]
# We can implement this using matrix multiplication carefully.
in_flow = (W.T @ p) # total inflow to each state i = sum_j W[j,i]*p_j
out_flow = (W @ p) # outflow from each state i = sum_j W[i,j]*p_i
master_term = in_flow - out_flow
# Feedback control u_i(t):
# Drive p(t) toward uniform distribution [1/3, 1/3, 1/3]
u = alpha * (1/3 - p)
# Sum up final ODE derivative
dpdt = master_term + u
return dpdt
def compute_entropy(p):
“”"
Compute Shannon entropy H(p) = - sum_i p_i log p_i
Safely handle p_i=0 by ignoring those terms (or adding a small epsilon).
“”"
p = np.asarray(p)
# For numerical stability, we can filter out p_i ~ 0 to avoid log(0).
mask = (p > 1e-12)
return -np.sum(p[mask] * np.log(p[mask]))
def run_dem_simulation(
p0=None,
W=None,
alpha=5.0,
t_span=(0, 10),
num_points=200
):
“”"
Simulate the DEM ODE with feedback for 3 states, then return time,
probabilities, and entropies.
Parameters
----------
p0 : array_like of shape (3,), optional
Initial distribution. If None, defaults to [0.8, 0.1, 0.1].
W : 2D array of shape (3, 3), optional
Transition-rate matrix. If None, a simple example is used.
alpha : float
Feedback gain for controlling the distribution.
t_span : (float, float)
Start and end time for ODE integration.
num_points : int
Number of time points to record for output.
Returns
-------
t_eval : np.ndarray
Time points at which solution is recorded.
p_sol : np.ndarray of shape (3, len(t_eval))
Probability distribution at each time point.
H_vals : np.ndarray
Entropy values at each time point.
"""
# Default initial distribution
if p0 is None:
p0 = np.array([0.8, 0.1, 0.1])
# Default transition-rate matrix
# W[i,j] is the rate from i -> j;
# Let's define something simple with no explicit time dependence:
if W is None:
# For example:
# - state 1 transitions to state 2 with rate 1.0, to state 3 with rate 0.5
# - state 2 transitions to state 1 with rate 0.3, to state 3 with rate 0.4
# - state 3 transitions to state 1 with rate 0.2, to state 2 with rate 0.1
# Diagonal entries can be set so that row sums do not necessarily matter here,
# because we handle inflow and outflow explicitly in the ODE.
W = np.array([
[0.0, 1.0, 0.5], # from state 1 to {1,2,3}
[0.3, 0.0, 0.4], # from state 2 to {1,2,3}
[0.2, 0.1, 0.0] # from state 3 to {1,2,3}
])
# Time grid for evaluating solution
t_eval = np.linspace(t_span[0], t_span[1], num_points)
# ODE solver call
sol = solve_ivp(
fun=lambda t, p: dynamic_entropy_model(t, p, W, alpha),
t_span=t_span,
y0=p0,
t_eval=t_eval
)
p_sol = sol.y # shape = (3, num_points)
# Compute entropy at each time point
H_vals = np.array([compute_entropy(p_sol[:, i]) for i in range(num_points)])
return sol.t, p_sol, H_vals
def plot_results(t_eval, p_sol, H_vals):
“”"
Generate plots: (1) p_i(t) vs time, (2) H(t) vs time.
“”"
fig, axs = plt.subplots(2, 1, figsize=(8, 6), sharex=True)
# Plot probabilities
axs[0].plot(t_eval, p_sol[0, :], label='p1(t)')
axs[0].plot(t_eval, p_sol[1, :], label='p2(t)')
axs[0].plot(t_eval, p_sol[2, :], label='p3(t)')
axs[0].set_ylabel('Probability')
axs[0].set_title('State Probabilities Over Time')
axs[0].legend(loc='best')
axs[0].grid(True)
# Plot entropy
axs[1].plot(t_eval, H_vals, 'r-', label='H(t)')
axs[1].set_xlabel('Time')
axs[1].set_ylabel('Entropy')
axs[1].set_title('Time-Dependent Entropy')
axs[1].legend(loc='best')
axs[1].grid(True)
plt.tight_layout()
plt.show()
def main():
# Run simulation
t_eval, p_sol, H_vals = run_dem_simulation(
p0=[0.8, 0.1, 0.1], # initial distribution
alpha=5.0, # feedback gain
t_span=(0, 10) # simulate from t=0 to t=10
)
# Plot results
plot_results(t_eval, p_sol, H_vals)
if name == “main”:
main()
Theoretical Proof Considerations
The numerical results strongly suggest that the system converges to uniform probability. However, to rigorously prove this:
- Fixed-Point Analysis: Show that the equilibrium distribution satisfies p_i^* = 1/3 for all i.
- Stability Analysis: Compute the Jacobian matrix and verify that its eigenvalues indicate global stability.
- Lyapunov Function Approach:
- Define entropy as a Lyapunov function:
[
V(p) = -\sum_i p_i \ln p_i.
] - Show that its time derivative \frac{dH}{dt} is non-negative, ensuring monotonic increase in entropy.
- Define entropy as a Lyapunov function:
This would constitute a formal proof that entropy increases under this control law and that the system stabilizes at the uniform distribution.
Extensions & Future Work
- Time-Varying W(t): Model changing external conditions (e.g., seasonal effects in biological systems).
- Different Feedback Laws: Instead of pushing to uniform, target a specific \mathbf{p}^*.
- Thermodynamic Cost Analysis: Track energy or entropy reservoirs.
- Higher-Dimensional Extensions: Generalize to N states.
- AI & RL Integration: Use reinforcement learning to optimize feedback control dynamically.
Conclusion
This Python framework encapsulates the essence of a Dynamic Entropy Model, demonstrating how:
- A small discrete-state system evolves probabilistically.
- Feedback control shapes probability distributions.
- Entropy changes over time, illustrating control effects.
This framework could be extended to more complex adaptive systems in physics, economics, and artificial intelligence.
DEM Simulation Analysis
Initial Probabilities | Final Probabilities | Initial Entropy | |
---|---|---|---|
State 1 | 0.8 | 0.262566406921853 | 0.639031859650177 |
State 2 | 0.1 | 0.34801296237766205 | 0.639031859650177 |
State 3 | 0.1 | 0.3699716553523358 | 0.639031859650177 |
The analysis of the simulation results shows:
- Probability Evolution:
- Initially, the system starts highly skewed ([0.8,0.1,0.1][0.8,0.1,0.1]).
- Over time, the feedback control shifts the probabilities closer to an even distribution ([0.26,0.35,0.37][0.26,0.35,0.37]), though not exactly uniform.
- Entropy Changes:
- The initial entropy was 0.639, which is relatively low due to the uneven distribution.
- The final entropy increased to 1.086, showing that the system evolved toward a higher-disorder state.
- The entropy increased by 0.447, confirming that the control mechanism drives the system toward a more uniform, higher-entropy state.
Visual Observations
- The probability curves show a smooth transition, with each state gradually approaching a balanced level.
- The entropy curve increases over time, matching theoretical expectations.
Conclusion
The results support the hypothesis that the feedback controller increases entropy and stabilizes the system near an even probability distribution. However, the final state is not perfectly uniform, indicating that further tuning (e.g., a higher ααvalue) might improve convergence.