Reaction Diffusion Results

The BZ reaction is implemented in a reaction diffusion framework as dicussed in the background. The system we are solving is then

Parameter values are as listed in the background. The diffusion constants are D_x=D_y=5e-3, D_z=1e-4.

Numerical method 1

The first implementation of the reaction diffusion system uses finite differences to discretize the spatial laplacian. Time stepping is implemented using the built-in Matlab stiff ODE solver ode15s. The primary difficulty in numerically solving these equations is the very short time scale introduced in the peaks of the concentrations relative to the time scale of the relaxation phase.

The equations are solved on a periodic square domain, [-π, π]x[-π, π], with N=48 points in each direction.

The initial conditions are similar to the Brusselator intitial conditions. Two points out of phase are selected and used as initial conditions for a bulge in the center of the domain and the remainder of the domain.

The result is shown in the embedded video below. The code can be found in https://git.uwaterloo.ca/tghill/chem-oscillations/tree/master/code/fkn_solver_finite_difference.

The behaviour of the system can be compared to the ODE model. Tracking the concentrations at the edge of the domain and at the center of the domain, the behaviour is overall similar to the ODE model.

As expected, the concentrations at the center and the edge are out of phase from each other. The behaviour is qualitatively the same as the ODE model. If the model equations are integrated longer, the phase lag reduces over time due to diffusion.

Clearly we want to do better than this solution.

Numerical method 2

A better implementation is achieved by taking the same approach as for the Brusselator model. Modifying the Brusselator code to work for a three component system, we arrive at the code found in https://git.uwaterloo.ca/tghill/chem-oscillations/tree/master/code/fkn_solver_spectral.

Unlike the Brusselator system, the FKN system is extremely stiff. Since we now use our own time stepping method, we must be careful with the time step. From trial and error, the time step must be on the order of dt=2e-7 for stability. Larger, and the system becomes unstable leading to nan concentrations. However, each cycle takes on the order of 5 units of time to complete. So for even the shortest useful simulation, we require approximately 10,000,000 time steps. This strongly limits the feasible spatial resolution.

The equations are solved on the same domain, with N=256 points in each dimension, for 4 units of time. The modelled X, Y, and Z concentrations are shown in the following videos. Parameters are taken as in the ODE model, with f=0.5. The diffusion constants are ru=rv=rw=0.01. This higher diffusion smooths out numerical artifacts. Smaller diffusion constants lead to NaN or unphysical solutions.

There are some artifacts in the solution, particular in the X concentration since it is the most rapidly varying, although the timeseries of the concentrations (taken at the center of the domain) matches well with the ODE results.

How to improve

With the number of required time steps, solving these equations over a long time interval at high spatial resolution reasonably requires parallelization. The reaction component can be easily decoupled and run in parallel, which would significantly speed up the code. With this model, it is worthwhile to spend time thinking about the implementation.