SimpleImplicitTauLeaping solver

[sivasathyaseeelan]

using JumpProcesses, DiffEqBase
using Test, LinearAlgebra, Statistics
using StableRNGs, Plots
rng = StableRNG(12345)

# Example system from Cao et al. (2007)
function define_stiff_system()
    # Reactions: S1 -> S2, S2 -> S1, S2 -> S3
    # Rate constants
    c = (1000.0, 1000.0, 1.0)
    
    # Define MassActionJump
    # Reaction 1: S1 -> S2
    reactant_stoich1 = [Pair(1, 1)]  # S1 consumed
    net_stoich1 = [Pair(1, -1), Pair(2, 1)]  # S1 -1, S2 +1
    # Reaction 2: S2 -> S1
    reactant_stoich2 = [Pair(2, 1)]  # S2 consumed
    net_stoich2 = [Pair(1, 1), Pair(2, -1)]  # S1 +1, S2 -1
    # Reaction 3: S2 -> S3
    reactant_stoich3 = [Pair(2, 1)]  # S2 consumed
    net_stoich3 = [Pair(2, -1), Pair(3, 1)]  # S2 -1, S3 +1
    
    jumps = MassActionJump([c[1], c[2], c[3]], [reactant_stoich1, reactant_stoich2, reactant_stoich3], 
                          [net_stoich1, net_stoich2, net_stoich3])
    
    return jumps
end

# Simulation parameters
u0 = [100, 0, 0]  # Initial conditions: S1=100, S2=0, S3=0
tspan = (0.0, 5.0)
p = nothing  # No additional parameters needed

# Define problem
prob = DiscreteProblem(u0, tspan, p)
jump_prob = JumpProblem(prob, PureLeaping(), define_stiff_system())

# Solve with SimpleImplicitTauLeaping
alg = SimpleImplicitTauLeaping()
sol = solve(jump_prob, alg; saveat=0.01)
plot(sol)

using JumpProcesses, DiffEqBase
using Test, LinearAlgebra
using StableRNGs, Plots
rng = StableRNG(12345)

# Reversible isomerization system from Cao et al. (2004)
function define_reversible_isomerization()
    c = (1000.0, 1000.0)
    reactant_stoich1 = [Pair(1, 1)]
    net_stoich1 = [Pair(1, -1), Pair(2, 1)]
    reactant_stoich2 = [Pair(2, 1)]
    net_stoich2 = [Pair(1, 1), Pair(2, -1)]
    jumps = MassActionJump([c[1], c[2]], [reactant_stoich1, reactant_stoich2], 
                          [net_stoich1, net_stoich2])
    return jumps
end

# Simulation parameters
u0 = [1000, 0]
tspan = (0.0, 0.1)
p = nothing
prob = DiscreteProblem(u0, tspan, p)
jump_prob = JumpProblem(prob, PureLeaping(), define_reversible_isomerization())

# Solve with both solvers
alg_newton = SimpleImplicitTauLeaping(epsilon=0.05, solver=NewtonImplicitSolver())
sol_newton = solve(jump_prob, alg_newton)
plot(sol_newton)

alg_trapezoidal = SimpleImplicitTauLeaping(epsilon=0.05, solver=TrapezoidalImplicitSolver())
sol_trapezoidal = solve(jump_prob, alg_trapezoidal)
plot(sol_trapezoidal)

NewtonImplicitSolver

TrapezoidalImplicitSolver

Checklist

• Appropriate tests were added
• Any code changes were done in a way that does not break public API
• All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC.
• Any new documentation only uses public API

Additional context

Add any other context about the problem here.