Solving Differential Equations

Available packages for ODEs and PDEs

The DifferentialEquations package is a very powerful toolkit for solving differential equations. It has multiple, state of the art solvers for ordinary differential equations. It can also handle stochastic differential equations, random differential equations, delay differential equations, differential algebraic equations, stochastic differential algebraic equations etc.

Should you wish to solve partial differential equations, you will need to convert the PDE to a set of ODEs. This can be done via the method of lines (MethodOfLines.jl), or via several other advanced tools in Julia, including Trixi.jl or several other finite element methods see JuliaFEM.

Here we shall focus on the most common applications in process modelling, solving systems of ordinary differential equations. In a later section, we shall look at using the method of lines to solve partial differential equations, along with specifying the equations in symbolic form.

System of ODEs - Solving a Plug Flow Reactor Model

For our ODE example, we take a simple PFR model and a kinetic model.

\[ A + B \rightarrow C\]

with the reaction proceeding in the gas phase over a catalyst with the following rate equation:

\[-r_A = k_o exp \Big[\frac{E_a}{R}(1/T - 1/T_o)\Big] P_A P_B\] \[[-r_A] = kmol/h/kg\]

We also have the following data:

Parameter	Value
\(k_o\)	\(0.5 \frac{kmol}{h.kg.bar^2}\)
\(E_a\)	\(100\ J/mol\)
\(T_o\)	\(500\ K\)
\(R\)	\(8.314\ J/mol/K\)

The feed stream entering the reactor has 60kmol/h of A and 40kmol/h of B. The reactor is isothermal at 510K. The inlet pressure is 10bar(abs) and the pressure drop is 1bar. Calculate the composition profile for the reactor.

In order to solve the problem, we need to tell the solver what the set of ODEs is for each of the molar flows. We already know the equation for the rate of change for A as it flows over the catalyst. Stoichiometry gives us the rate of change for B and C.

We package this information into a function, rates(u, p, t), which is what we will supply to the solver.

The differential equation solver will accept two forms of functions:

rates(u, p, t) is a function that takes the current values of the dependent variable(s), model parameter(s), p and the current value of the independent variable, t, and then returns the differential(s) of the dependent variable(s) at that point.

u can be a vector or a scalar. p should be any iterable type, but is usually passed as a tuple, since this will give faster code than passing a vector.

The alternative form, called the in-place form, is rates!(du, u, p, t). In this case, du is a scalar or vector that holds the derivative(s) and is updated in the function. This in-place version typically results in faster code, as there are fewer memory allocation needed. This should be your default form for specifying ordinary differential equations. The exclamation mark (bang, for our American friends) is just a Julia convention to indicate that this is a mutating function, i.e. it modifies one or more of its parameters. Again by convention, this is will usually be the first parameter, du. You do not have to add the exclamation mark to the function name, but this is definitely the preferred thing to do.

using DifferentialEquations, Plots

"""
The reactions rates for A + B -> C
    u:  array with the molar flows of A, B and C
    p:  tuple with P and T
    t:  position in the reactor 0..L
"""
function rates!(du, u, p, t)
    k₀ = 0.50
    Eₐ = 100.0
    T₀ = 500.0
    R = 8.314
    
    # Convert molar flows to molar fractions
    y = u ./ sum(u)

    # Unpack the tuple to get pressure and temperature
    P₀, ΔP, T, tmax = p

    P = P₀ - t/tmax*ΔP
    k = k₀ * exp(-Eₐ/R * (1/T - 1/T₀))
    rA = -k * y[1] * y[2] * P^2
    rB = rA
    rC = -rA

    du[1] = rA
    du[2] = rB
    du[3] = rC
end

# Set up the input values and pack into a tuple of parameters
T = 510.0
P = 10.0
ΔP = 1.0
tmax = 20.0
p = (P, ΔP, T, tmax)

# The starting values for the initial value problem
u₀ = [60.0, 40.0, 0.0]

# The span over which to integrate - here the mass of catalyst
tspan = (0.0, tmax)

# Create the ODEProblem object
prob = ODEProblem(rates!, u₀, tspan, p)

# Solve with the default solver
sol = solve(prob)

# Plot the result - the molar flows vs the mass of catalyst
plot(sol, leg=:right)
savefig("img/ODEmol.svg")

The plot of sol.u vs sol.t, generated by `plot(sol)

The default plot recipe (defined by the DifferentialEquations package) will plot each of the dependent variables against the independent variable. All we need to do is call plot(sol). Here, the default placement of the legend was interfering with the lines, so we move it to the right centre, by adding leg = :right) to the call to plot().

Now we have solved the set of three differential equations and we have the molar flows of each component at various points along the length of the reactor. The points were set by the variable step-size solver. There are many options you can pass the solver, including fixed steps at which to record the values. See the documentation for details.

Part of what is returned in the sol structure is an interpolator. If we want to have the molar flows at any point in the reactor, we need only call sol(x) to get the values of the dependent variables (sol.u) at a value x of the independent variable (sol.t). sol is a callable structure, sometimes called a functor.

sol(5.123)
# 3-element Vector{Float64}:
#  28.55013685340714
#   8.550136853407144
#  31.44986314659286

The original question however was for the composition, which is generally interpreted as the fractions of the components, not their flows. Let’s convert the results to molar fractions:

# Now, we want to plot the molar fractions...
# Copy over the independent values (sol.t) into our x array (optional)
x = sol.t

# For each entry in sol.u, sum up the molar flows - we get an array of values
# with the total molar flow at every point in the reactor
totflows = sum.(sol.u)

# Now divide the individual flows with the total. 
fracs = sol.u ./ totflows

# We get the same structure as sol.u - an array of arrays.
# We need to flatten this to plot it.
y = hcat(fracs...)'
plot(x,y, labels=["A" "B" "C"])

The plot of molar fractions in the reactor profile

The last bit of code in this part will need some explanation. As mentioned above, the solution returned by the ODE solver is a structure, sol with multiple fields, including sol.t for the independent variable. This will hold the values for each step taken by the solver. The results for the dependent variables, u[...] are in sol.u, for each the steps in t. This means that sol.u is an array of arrays:

sol.u
# 17-element Vector{Vector{Float64}}:
#  [60.0, 40.0, 0.0]
#  [59.99858583839893, 39.99858583839893, 0.001414161601070281]  
#  [59.984446697061045, 39.984446697061045, 0.015553302938956699]
#  [59.84330258351507, 39.84330258351507, 0.15669741648493382]   
#  [58.45642225894246, 38.45642225894246, 1.543577741057545]     
#  ⋮
#  [21.04774569180567, 1.0477456918056731, 38.952254308194334]
#  [20.551956292857035, 0.5519562928570385, 39.44804370714297]
#  [20.282897789407816, 0.2828977894078193, 39.71710221059219]
#  [20.20298976079487, 0.20298976079487502, 39.79701023920513]

We can to convert these to molar fractions, which means we need to divide each entry of the sub-arrays by the sum of that sub-array. This is very easy in Julia. We just do this:

totflows = sum.(sol.u)
fracs = sol.u ./ totflows

We broadcast the sum() function over sol.u by inserting a period between the function name and the first parenthesis - sum.(). This means the function is applied individually to each element in sol.u. Since each entry of sol.u is an array, we therefore get the an array of the total molar flows at each point in the reactor.

Now we want to divide each entry in each sub-array in sol.u with this sum to get the molar fractions at each point. We do this by broadcasting again. For functions the period is inserted between the function name and the parenthesis, but for operators, the period is put in front of the operator. The broadcast, or element-wise, division will divide each entry in fracs (an array of molar flows) by the corresponding entry in totflows. The result in an array of arrays, just like sol.u, but with the molar fractions instead of molar flows.

totflows = sum.(sol.u)
# 17-element Vector{Float64}:
#  100.0
#   99.99858583839894
#   99.98444669706105
#   99.84330258351507
#   98.45642225894247
#    ⋮
#   61.04774569180567
#   60.55195629285704
#   60.282897789407826
#   60.20298976079488

fracs = sol.u ./ totflows
# 17-element Vector{Vector{Float64}}:
#  [0.6, 0.4, 0.0]
#  [0.5999943432736005, 0.3999915149104007, 1.4141815998834358e-5]
#  [0.5999377771105296, 0.3999066656657945, 0.0001555572236758087]
#  [0.5993722266294071, 0.3990583399441108, 0.0015694334264820865]
#  [0.5937288895710717, 0.39059333435660765, 0.015677776072320636]
#  ⋮
#  [0.34477515022525834, 0.01716272533788762, 0.6380621244368541]
#  [0.339410277571518, 0.009115416357277122, 0.6514743060712049]
#  [0.33646189107006863, 0.0046928366051030585, 0.6588452723248283]
#  [0.3355811703217333, 0.0033717554826000534, 0.6610470741956666]

The plotting library expects you to pass the x-values as an array and when there are multiple y-values, these should be in a matrix where each column is a series. We therefore need to convert the array of arrays into a matrix.

This is done with y = hcat(fracs...)', which uses two new Julia concepts.

The first is splatting. If a function, say f(a, b, c) takes three scalar parameters, we can splat an array into the function by calling f(d...), where d is an array. The first entry in d will be passed into a, the second into b and the rest of the array into c. If d had only three entries, c will be a scalar, else c will be a tuple.

The opposite of splatting is slurping, when multiple scalars are slurped up into a single array parameter. It uses the same syntax of an ellipsis.

What happens here is that the hcat() function (horizontal concatenation) expects a series of values that will be concatenated as the columns of a new array. We pass it the fracs variable, which is an array of arrays, and we splat it. This means that hcat() gets a series of inputs, one for each row in fracs. These are then concatenated as follows:

hcat(fracs...)
# 3×17 Matrix{Float64}:
#  0.6  0.599994    0.599938     0.599372    0.593729   0.574037   0.540935  0.502011  …  0.391779   0.369135  0.354054   0.344775   0.33941     0.336462    0.335581
#  0.4  0.399992    0.399907     0.399058    0.390593   0.361056   0.311402  0.253017     0.0876678  0.053702  0.0310817  0.0171627  0.00911542  0.00469284  0.00337176       
#  0.0  1.41418e-5  0.000155557  0.00156943  0.0156778  0.0649074  0.147664  0.244972     0.520554   0.577163  0.614864   0.638062   0.651474    0.658845    0.661047

Why use hcat() and not vcat()? Remember that a 1D array in Julia is considered to be a column vector. This means fracs is a column vector where each entry is also a column vector. If we vcat fracs, these individual column vectors get appended to each other in one long column vector and we lose the distinction between the individual points along the reactor:

vcat(fracs...)
# 51-element Vector{Float64}:
#  0.6
#  0.4
#  0.0
#  0.5999943432702689
#  0.3999915149054033
#  1.4141824327869601e-5
#  0.5999377767074147
#  ⋮
#  0.33587562016618194
#  0.00381343024927297
#  0.6603109495845451
#  0.33464693236081094
#  0.0019703985412164293
#  0.6633826690979727

But the hcat(fracs...) call horizontally concatenated the sub-arrays (column vectors), making each series (species) a row of the matrix. We need each series (species) to be a column. We get that by transposing the result with the ' appended to the end. A' is the transposed version of a real matrix, A

Important

A' is technically the adjoint matrix. This means it is the transpose of the matrix of complex conjugates. For real values, however, this makes no difference. For complex values, you need to call transpose(), rather than use the adjoint operator, if you want the transpose.

y = hcat(fracs...)'
# 17×3 adjoint(::Matrix{Float64}) with eltype Float64:
#  0.6       0.4         0.0
#  0.599994  0.399992    1.41418e-5
#  0.599938  0.399907    0.000155557
#  0.599372  0.399058    0.00156943
#  0.593729  0.390593    0.0156778
#  ⋮
#  0.344775  0.0171627   0.638062
#  0.33941   0.00911542  0.651474
#  0.336462  0.00469284  0.658845
#  0.335581  0.00337176  0.661047

Finally we get:

y = hcat(fracs...)'
# 17×3 adjoint(::Matrix{Float64}) with eltype Float64:
#  0.6       0.4         0.0
#  0.599994  0.399992    1.41418e-5
#  0.599938  0.399907    0.000155558
#  0.599372  0.399058    0.00156954
#  0.593725  0.390587    0.015688
#  0.573817  0.360726    0.0654572
#  0.540162  0.310243    0.149596
#  ⋮
#  0.366919  0.0503788   0.582702
#  0.352276  0.0284133   0.619311
#  0.343461  0.0151917   0.641347
#  0.338506  0.00775872  0.653735
#  0.335876  0.00381343  0.660311
#  0.334647  0.0019704   0.663383

Note

Note that the type returned is an adjoint(::Matrix{Float64}). This is a lazy structure. It just maps onto the original matrix, rather than recalculate it, for better efficiency.

Now we have each series as a column and we pass the series labels as a row vector (with spaces between the entries, not commas), since each series must also be a column in the labels parameter.