You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Copy file name to clipboardExpand all lines: docs/src/model_creation/examples/smoluchowski_coagulation_equation.md
+66-50Lines changed: 66 additions & 50 deletions
Display the source diff
Display the rich diff
Original file line number
Diff line number
Diff line change
@@ -4,72 +4,89 @@ This tutorial shows how to programmatically construct a [`ReactionSystem`](@ref)
4
4
The Smoluchowski coagulation equation describes a system of reactions in which monomers may collide to form dimers, monomers and dimers may collide to form trimers, and so on. This models a variety of chemical/physical processes, including polymerization and flocculation.
5
5
6
6
We begin by importing some necessary packages.
7
-
```julia
7
+
```@example smcoag1
8
8
using ModelingToolkit, Catalyst, LinearAlgebra
9
-
usingDiffEqBase, JumpProcesses
9
+
using JumpProcesses
10
10
using Plots, SpecialFunctions
11
11
```
12
12
Suppose the maximum cluster size is `N`. We assume an initial concentration of monomers, `Nₒ`, and let `uₒ` denote the initial number of monomers in the system. We have `nr` total reactions, and label by `V` the bulk volume of the system (which plays an important role in the calculation of rate laws since we have bimolecular reactions). Our basic parameters are then
13
-
```julia
14
-
## Parameter
15
-
N =10# maximum cluster size
16
-
Vₒ = (4π/3)*(10e-06*100)^3# volume of a monomers in cm³
The [Smoluchowski coagulation equation](https://en.wikipedia.org/wiki/Smoluchowski_coagulation_equation) Wikipedia page illustrates the set of possible reactions that can occur. We can easily enumerate the `pair`s of multimer reactants that can combine when allowing a maximal cluster size of `N` monomers. We initialize the volumes of the reactant multimers as `volᵢ` and `volⱼ`
26
35
27
-
```julia
36
+
```@example smcoag1
28
37
# possible pairs of reactant multimers
29
38
pair = []
30
39
for i = 2:N
31
-
push!(pair, [1:integ(i/2) i .- (1:integ(i/2))])
40
+
halfi = floor(Int, i/2)
41
+
push!(pair, [(1:halfi) (i .- (1:halfi))])
32
42
end
33
43
pair = vcat(pair...)
34
-
vᵢ =@view pair[:,1] # Reactant 1 indices
35
-
vⱼ =@view pair[:,2] # Reactant 2 indices
36
-
volᵢ = Vₒ*vᵢ # cm⁻³
37
-
volⱼ = Vₒ*vⱼ # cm⁻³
44
+
vᵢ = @view pair[:, 1] # Reactant 1 indices
45
+
vⱼ = @view pair[:, 2] # Reactant 2 indices
46
+
volᵢ = Vₒ * vᵢ # cm⁻³
47
+
volⱼ = Vₒ * vⱼ # cm⁻³
38
48
sum_vᵢvⱼ = @. vᵢ + vⱼ # Product index
49
+
nothing #hide
39
50
```
40
51
We next specify the rates (i.e. kernel) at which reactants collide to form products. For simplicity, we allow a user-selected additive kernel or constant kernel. The constants(`B` and `C`) are adopted from Scott's paper [2](https://journals.ametsoc.org/view/journals/atsc/25/1/1520-0469_1968_025_0054_asocdc_2_0_co_2.xml)
41
-
```julia
52
+
```@example smcoag1
42
53
# set i to 1 for additive kernel, 2 for constant
43
54
i = 1
44
-
if i==1
45
-
B =1.53e03# s⁻¹
46
-
kv =@. B*(volᵢ + volⱼ)/V # dividing by volume as its a bi-molecular reaction chain
55
+
if i == 1
56
+
B = 1.53e03 # s⁻¹
57
+
58
+
# dividing by volume as it is a bimolecular reaction chain
59
+
kv = @. B * (volᵢ + volⱼ) / V
47
60
elseif i==2
48
-
C =1.84e-04# cm³ s⁻¹
49
-
kv =fill(C/V, nr)
61
+
C = 1.84e-04 # cm³ s⁻¹
62
+
kv = fill(C / V, nr)
50
63
end
64
+
nothing #hide
51
65
```
52
-
We'll store the reaction rates in `pars` as `Pair`s, and set the initial condition that only monomers are present at ``t=0`` in `u₀map`.
53
-
```julia
54
-
# unknown variables are X, pars stores rate parameters for each rx
66
+
We'll set the parameters and the initial condition that only monomers are present at ``t=0`` in `u₀map`.
67
+
```@example smcoag1
68
+
# k is a vector of the parameters, with values given by the vector kv
69
+
@parameters k[1:nr] = kv
70
+
71
+
# create the vector of species X_1,...,X_N
55
72
t = default_t()
56
-
@species k[1:nr] (X(t))[1:N]
57
-
pars =Pair.(collect(k), kv)
73
+
@species (X(t))[1:N]
58
74
59
75
# time-span
60
76
if i == 1
61
-
tspan = (0. ,2000.)
77
+
tspan = (0.0, 2000.0)
62
78
elseif i == 2
63
-
tspan = (0. ,350.)
79
+
tspan = (0.0, 350.0)
64
80
end
65
81
66
82
# initial condition of monomers
67
83
u₀ = zeros(Int64, N)
68
84
u₀[1] = uₒ
69
-
u₀map =Pair.(collect(X), u₀) # map variable to its initial value
85
+
u₀map = Pair.(collect(X), u₀) # map species to its initial value
86
+
nothing #hide
70
87
```
71
88
Here we generate the reactions programmatically. We systematically create Catalyst `Reaction`s for each possible reaction shown in the figure on [Wikipedia](https://en.wikipedia.org/wiki/Smoluchowski_coagulation_equation). When `vᵢ[n] == vⱼ[n]`, we set the stoichiometric coefficient of the reactant multimer to two.
We now convert the [`ReactionSystem`](@ref) into a `ModelingToolkit.JumpSystem`, and solve it using Gillespie's direct method. For details on other possible solvers (SSAs), see the [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/types/jump_types/) documentation
0 commit comments