The Matrix Market

Home Search Browse Resources

MVMBWM: Brusselator Wave Model in Chemical Reaction

from the NEP Collection

Matrix Generator MVMBWM
Source: Y. Saad, University of Minnesota
Discipline: Chemical engineering
Language: Fortran
Output format: matrix-vector multiply

This problem models the concentration waves for reaction and transport interaction of chemical solutions in a tubular reactor. The concentrations x(t,z) and y(t,z) of two reacting and diffusing components are modeled by the system

dx/dt=delta_1/L^2 d^2x/dz^2+f(x,y); dy/dt = delta_2/L^2 d^2y/dz^2 +g(x,y)

with the initial conditions x(0,z)=x_0(z); y(0,z)=y_0(z), and the Dirichlet boundary conditions x(t,0)=x(t,1)=x^*; y(t,0)=y(t,1)=y^*, where 0 <= z <= 1 is the space coordinate along the tube, and t is time. Raschman et al. considered in particular the so-called Brusselator wave model in which

f(x,y)=alpha-(beta+1)*x +x^2*y; g(x,y)=beta*x-x^2*y

. Then, the above system admits the trivial stationary solution x^*=alpha; y^*=beta/alpha. In this problem one is primarily interested in the existence of stable periodic solutions to the system as the bifurcation parameter L varies. This occurs when the eigenvalues of largest real parts of the Jacobian of the right hand side of (1) and (2) are evaluated at the steady station solution, is purely imaginary. For the purpose of verifying this fact numerically, one first needs to discretize the equations with respect to the variable z and compute the eigenvalues with largest real parts of the resulting discrete Jacobian.

If we discretize the interval [0,1] using n interior points with the mesh size h=1/(m+1), then the discretized Jacobian of the system is a 2x2 block matrix of the form


where T = tridiag{1,-2,1}, tau_1=1/h^2 delta_1/L^2 and tau_2=1/h^2 delta_2/L^2.

The exact eigenvalues of J are known since there exists a quadratic relation between the eigenvalues of the matrix A and those of the classical difference matrix T = tridiag{1,-2,1}. The following segment is the Matlab M-file for computing the 2m eigenvalues of A: mvmbwm M-file segment.

Figure 1 shows the 32 rightmost eigenvalue distribution of 200 by 200 of the matrix A (m=100) corresponding to the set of parameters

delta_1=0.008; delta_2=delta_1/2=0.004;  alpha=2; beta=5.45; L=0.51302

as used in Saad's book.

The following is the FORTRAN calling sequence for forming matrix-vector AX or ATX: Fortran calling sequence. In addition, the parameters delta_1, delta_2, L, alpha and beta are set as internal parameters.


Norder of the matrix (2m)
ALPHAconstant in reaction term for x
BETAcoefficient in reaction term for y
DELT1diffusion coefficient for x
DELT2diffusion coefficient for y
Lbifurcation parameter (L2 divides DELT1 and DELT2)


Matrix Instances:

The Matrix Market is a service of the Mathematical and Computational Sciences Division / Information Technology Laboratory / National Institute of Standards and Technology.

[ Home ] [ Search ] [ Browse ] [ Resources ]

Last change in this page: Wed Sep 22 13:37:30 US/Eastern 2004 [Comments: ]