ReactionDiffusionEqn is a Julia package that can find numerical solutions to one-dimensional reaction diffusion equations based on the following input:
- Diffusion Coefficient (D)
- Reaction Term (f(x, t, u))
- Initial Condition (u₀)
- Boundary Conditions (Dirichlet, Neumann, Convective Heat (aka Robin), or Periodic)
bc = Dirichlet(ū) # boundary condition must be specified
bc = Neumann(∂ū) # boundary condition derivative must be specified
bc = Periodic() # nothing needs to be specified
bc = ConvectiveHeat(T̄₀, K̄) # ambient temperature (T̄₀) and thermal conductivity (K̄) must be specified- Number of Spatial Steps (Nₓ - number of spatial discretization points)
- Space Time (space-time over which solution to PDE is approximated)
st = SpaceTime(L, tf) # L is the spatial extent and tf is the time span of simulation)- Sample Time (stores u every sample_time time steps)
The solution to the PDE is approximated numerically using a finite difference spatial and temporal discretization with the Crank-Nicolson Method:
Approximating this complex equation is as easy as passing eight (seven if periodic, see periodic example for further explanation) variables into the solve_rnx_diffn_eqn function. Below is an example with Dirichlet Boundary Conditions and the corresponding heat map:
D = 1.0
f(x::Float64, t::Float64, u::Float64) = 100 * (D * π^2 - 1.0) * (e ^ (-t) * sin(π * x))
u₀(x::Float64) = 100 * sin(π * x)
left_bc = Dirichlet(0.0)
right_bc = Dirichlet(0.0)
Nₓ = 100
st = SpaceTime(1.0, 1.0)
sample_time = 0.1
t, x, u = solve_rxn_diffn_eqn(left_bc, right_bc, f, u₀, D, Nₓ, st, sample_time)
draw_heat_map(t, x, u)
ReactionDiffusionEqn can also handle mixed boundary conditions, such as Dirichlet-Neumann, Dirichlet-ConvectiveHeat, Neumann-ConvectiveHeat, etc.
The following functions are also available to aid in visualizations:
draw_heat_map(t, x, u)As shown in the example above, the draw_heat_map function outputs a heat map using PyPlot.
gif_maker(t, x, u, st, sample_time)The gif_maker function plots the approximated u at each timestep.
Use the ReactionDiffusionEqn package to determine how the temperature will change over time in an insulated, one-dimensional rod when it is heated at a particular location while the ends of the rod are held at a fixed temperature.
Imagine there is a one dimensional lake with a town at one end and a forest at the other as pictured above. The following equation can be used to model the fish population in the lake if people fish there:
The reaction term consists of the logistic growth model to simulate the population density of the fish and ch(x)p(u) to model the harvesting of the fish.
🐟 c is the harvesting rate.
🐟 h(x) is the harvesting distribution (i.e. how far from the town the person fishing drops a line in the lake).
🐟 p(u) is the probability of catching a fish based on the population density of the fish.
The changes in grass density can be modeled where a tree shades some grass and fauna eat the grass with the following equation to simulate an infinite field with Periodic boundary conditions:
The reaction term consists of the logistic growth model to simulate the density of the grass and - mb - hqb to model death of the grass.
🐘 g is the growth rate of the grass
🐘 m is the mortality rate of the grass
🐘 h is the biomass density of the fauna
🐘 q is the rate of consumption of the grass
When using Periodic boundary conditions, there is no need for a left and right boundary condition. Only one boundary condition is needed:
bc = Periodic()
t, x, u = solve_rxn_diffn_eqn(bc, f, u₀, D, Nₓ, st, sample_time)
The ReactionDiffusionEqn software can also handle mixed boundary conditions. In this example, Dirichlet and Convective Heat boundary conditions are mixed. Above is an insulated, one-dimensional rod, and this software can be used to evaluate the change in temperature as the rod is heated while one end of the rod is held at a fixed temperature.