The rest of this page has not been updated.

Version 1.5, January 2003

Selwyn Hollis
s.hollis@earthlink.net
http://www.math.armstrong.edu/faculty/hollis

Description

ReactionDiffusionLab is a Mathematica package that provides two functions, RDDensityPlots and RDSurfacePlots, for animating solutions of autonomous reaction-diffusion systems with two components,

or three components,

The domain is the unit square (, ), and boundary conditions are any combination of homogeneous Neumann type and (possibly nonhomogeneous) Dirichlet type.

Time stepping is done via an alternating direction implicit (ADI) method. (This is the primary improvement made to version 1.0.) Although the method is unconditionally stable for linear problems, a reasonably small time step is required for nonlinear problems.

At each step, solution values are stored in an matrix corresponding to a simple square grid with . Regardless of the boundary conditions, the first/last rows correspond to the top/bottom edges of the region, and the first/last entries in all rows correspond to the left/right edges of the region.

Limitations

The package works only on the problems described above. In particular there is no support for domains other than the unit square, spatially heterogeneous diffusion, Robin boundary conditions, or boundary conditions in which type depends on location. These are possible enhancements for future versions of the package.

Two-component Systems

The two functions are as follows in the two-component case:

In each of these, u0 and v0 are matrices containing initial values.

RDDensityPlots takes the options of ListDensityPlot, with the exception of FrameTicks. RDSurfacePlots takes the options of ListPlot3D. In addition, each has the following options:

PlotEvery	1
ReturnLast	False
BoundaryConditions	Neumann
BoundaryFunctions

The PlotRange option works in a nonstandard way: PlotRange -> {{a,b},{c,d}} provides separate vertical plot ranges for the two components.  With PlotRange -> {a,b} the same vertical plot range is used for both components. If no plot range is specified, or if PlotRange -> Automatic, then plot ranges are calculated from the minimum and maximum of the initial data. (With a nonnegative minimum these are {0,1.25 max}. ) This will rarely give an ideal choice, but it ensures that plot sequences always “stay still.”  A related matter: The default for the ClipFill option is None.

Examples

This loads the package (assuming the file is in AddOns : Applications : RDL) and clears anything that might cause a problem.

These are good values for the diffusion coefficients and .

Next define a matrix of initial values for each component.

Now sit back and watch. This is the Brusselator system with Neumann boundary conditions.

Click on the image to see a QuickTime movie. (140k)

Here’s the same thing in 3D:

Click on the image to see a QuickTime movie. (220k)

The same system but with Dirichlet boundary data for ... and in color.

Click on the image to see a QuickTime movie. (160k)

In this one, both components are held at their equilibrium values on the boundary.

Click on the image to see a QuickTime movie. (210k)

Three-component Systems

For the three-component case, just include a third entry wherever appropriate:

Examples

This is the vector field for an Oregonator system:

This time we’ll use a 25×25-point grid.

This defines initial values for each component.

The following uses the Dirichlet boundary condition and Neumann boundary conditions for and .

Click on the image to see a QuickTime movie. (92k)

Here’s the same thing in 3D:

Click on the image to see a QuickTime movie. (80k)

If you have the patience, run the following. It gives a nice animation of the oscillatory behavior of the solution.

Click on the image to see a QuickTime movie. (625k)