The Brusselator example: Continuation of a solution to a boundary value problem in a free parameter

Discretized solutions of PDEs can also be continued in CL_MATCONT. We illustrate this by continuing the equilibrium solution to a one-dimensional PDE. The curve type is called 'pde_1'. It is defined in the file pdf_1.m in the directory Equilibrium.

The Brusselator is a system of equations intended to model the Belusov-Zhabotinsky reaction. This is a system of reaction-diffusion equations that is known to exhibit oscillatory behavior. The unknowns are the concentrations $X(x,t),Y(x,t), A(x,t)$ and $B(x,t)$ of four reactants. Here $t$ denotes time and $x$ is a one-dimensional space variable normalized so that $x\in [0,1]$. The length $L$ of the reactor is a parameter of the problem. In our simplified setting $A$ and $B$ are constants.

The system is described by two partial differential equations:

\begin{displaymath}
\begin{array}{ccl}
\ensuremath{\frac{\partial X}{\partial...
...th{\frac{\partial^2Y}{\partial x^2}} + BX - X^2Y
\end{array}
\end{displaymath} (99)

with $x\in[0,1],\ t\geq 0$. Here $D_x,D_y$ are the diffusion coefficients of $X$ and $Y$. At the boundaries $x=0$ and $x=1$ Dirichlet conditions are imposed:
\begin{displaymath}
\left\{ \begin{array}{l}
X(0,t) = X(1,t) = A \\
Y(0,t) = Y(1,t) = \frac{B}{A}
\end{array}\right.
\end{displaymath} (100)

We are interested in equilibrium solutions $X(x)$ and $Y(x)$ to the system and their dependence on the parameter $L$.

The approximate equilibrium solution is:

\begin{displaymath}
\left\{ \begin{array}{ccl}
X(x) &=& A + 2\sin(\pi x) \\ 
...
... &=& \frac{B}{A} - \frac{1}{2}\sin(\pi x)
\end{array}\right.
\end{displaymath} (101)

The initial values of the parameters are: $A=2$, $B=4.6$, $D_x=0.0016$, $D_y=0.08$ and $L=0.06$. The initial solution (101) is not an equilibrium, but the continuer will try to converge to an equilibrium close to the initial solution. We use equidistant meshes. To avoid spurious solutions (solutions that are induced by the discretization but do not actually correspond to solutions of the undiscretized problem) one can vary the number of mesh points by setting the parameter $N$. If the same solution is found for several discretizations, then we can assume that they correspond to solutions of the continuous problem.

The second order space derivative is approximated using the well-known three-points difference formula: $\ensuremath{\frac{\partial^2f}{\partial x^2}} = \frac{1}{h^2}\left(f_{i-1} - 2f_i + f_{i+1}\right)$, where $h=\frac{1}{N+1}$, where $N$ is the number of grid points on which we discretize $X$ and $Y$. So $N$ is a parameter of the problem and $2N$ is the number of state variables (which is not fixed in this case).

The Jacobian is a sparse 5-band matrix. In the ode-file describing the problem the Jacobian is introduced as a sparse matrix. The Hessian is never computed as such but second order derivatives are computed by finite differences whenever needed. We note that MATLAB $6.5$ or $7$ does not provide sparse structures for 3-dimensional arrays.

bruss.m
1 function out = bruss
2 %
3 % Odefile of 1-d Brusselator model
4 %
5
6 out{1} = @init;
7 out{2} = @fun_eval;
8 out{3} = @jacobian;
9 out{4} = @jacobianp;
10 out{5} = [];%@hessians;
11 out{6} = [];%@hessiansp;
12 out{7} = [];
13 out{8} = [];
14 out{9} = [];
15
16
17
18 % --------------------------------------------------------------------------
19
20
21 function dfdt = fun_eval(t,y,N,L)
22
23 x = y(1:N);
24 y = y(N+1:2*N);
25
26 A = 2;
27 B = 4.6;
28 Dx = 0.0016;
29 Dy = 0.008;
30 x0 = A; x1 = A;
31 y0 = B/A; y1 = B/A;
32 L2 = L^2;
33 h = 1/(N+1);
34 cx = (Dx/L2)/(h*h);
35 cy = (Dy/L2)/(h*h);
36
37 dxdt = zeros(N,1);
38 dydt = zeros(N,1);
39
40 dxdt(1) = (x0-2*x(1)+x(2))*cx + A - (B+1)*x(1) + x(1)*x(1)*y(1);
41 dxdt(N) = (x(N-1)-2*x(N)+x1)*cx + A - (B+1)*x(N) + x(N)*x(N)*y(N);
42
43 dydt(1) = (y0-2*y(1)+y(2))*cy + B*x(1) - x(1)*x(1)*y(1);
44 dydt(N) = (y(N-1)-2*y(N)+y1)*cy + B*x(N) - x(N)*x(N)*y(N);
45
46 for i=2:N-1
47 dxdt(i) = (x(i-1)-2*x(i)+x(i+1))*cx + A - (B+1)*x(i) + x(i)*x(i)*y(i);
48 dydt(i) = (y(i-1)-2*y(i)+y(i+1))*cy + B*x(i) - x(i)*x(i)*y(i);
49 end
50
51 dfdt = [dxdt; dydt];
52
53 % --------------------------------------------------------------------------
54
55 function [tspan,y0,options] = init(N)
56 tspan = [0; 10];
57 A = 2;
58 B = 4.6;
59
60 y0 = zeros(2*N,1);
61
62 for i=1:N
63 y0(i) = A + 2*sin(pi*i/(N+1));
64 y0(N+i) = B/A - 0.5*sin(pi*i/(N+1));
65 end
66 handles = feval(@bruss);
67 options = odeset('Vectorized','on', 'Jacobian', handles(3), 'JacobianP', handles(4));
68
69 % --------------------------------------------------------------------------
70
71 function dfdxy = jacobian(t,y,N,L)
72 x = y(1:N);
73 y = y(N+1:2*N);
74 A = 2;
75 B = 4.6;
76 Dx = 0.0016;
77 Dy = 0.008;
78 x0 = A; x1 = A;
79 y0 = B/A; y1 = B/A;
80 L2 = L^2;
81 h = 1/(N+1);
82 cx = (Dx/L2)/(h*h);
83 cy = (Dy/L2)/(h*h);
84
85
86 %
87 % Sparse jacobian
88 %
89 A=zeros(2*N,3);
90 A(1:N-1,2)=cx;
91 A(1:N,3)=-2*cx -(B+1) + 2*x(1:N).*y(1:N);
92 A(1:N,4)=cx;
93
94 A(N+1:2*N,2) = cy;
95 A(N+1:2*N,3) = -2*cy -x(:).*x(:);
96 A(N+2:2*N,4) = cy;
97
98
99 A(1:N,1) = B - 2*x(:).*y(:);
100 A(N+1:2*N,5) = x(:).*x(:);
101
102 dfdxy = spdiags(A, [-N,-1:1,N] , 2*N, 2*N);
103 return
104
105 %
106 % Full matrix
107 %
108 dfxdx = zeros(N,N);
109 dfydy = zeros(N,N);
110 dfxdy = zeros(N,N);
111 dfydx = zeros(N,N);
112
113 for i=1:N
114 if i>1, dfxdx(i,i-1) = cx; end;
115 dfxdx(i,i) = -2*cx -(B+1) + 2*x(i)*y(i);
116 if i<N, dfxdx(i,i+1) = cx; end;
117 dfxdy(i,i) = x(i)*x(i);
118
119 if i>1, dfydy(i,i-1) = cy; end;
120 dfydy(i,i) = -2*cy -x(i)*x(i);
121 if i<N, dfydy(i,i+1) = cy; end;
122 dfydx(i,i) = B - 2*x(i)*y(i);
123 end
124
125 dfdxy = [ dfxdx, dfxdy; dfydx, dfydy ];
126
127 % --------------------------------------------------------------------------
128
129 function dfdp = jacobianp(t,y,N,L)
130 x = y(1:N);
131 y = y(N+1:2*N);
132 A = 2;
133 B = 4.6;
134 Dx = 0.0016;
135 Dy = 0.008;
136 x0 = A; x1 = A;
137 y0 = B/A; y1 = B/A;
138 L2 = L^2;
139 h = 1/(N+1);
140 cx = (Dx/L2)/(h*h);
141 cy = (Dy/L2)/(h*h);
142 kx = (-2/L)*cx;
143 ky = (-2/L)*cy;
144
145 Sx = zeros(N,1);
146 Sy = zeros(N,1);
147
148 Sx(1) = kx*(x0-2*x(1)+x(2));
149 Sy(1) = ky*(y0-2*y(1)+y(2));
150
151 Sx(N) = kx*(x(N-1)-2*x(N)+x1);
152 Sy(N) = ky*(y(N-1)-2*y(N)+y1);
153
154 i=2:N-1;
155 Sx(i) = kx*(x(i-1)-2*x(i)+x(i+1));
156 Sy(i) = ky*(y(i-1)-2*y(i)+y(i+1));
157
158 dfdp = [ zeros(2*N,1) [Sx;Sy] ];
159
160 % --------------------------------------------------------------------------
bruss.m
The model is implemented with 2 parameters: $N$ and $L$; the values of $A,B,D_x,D_y$ are hard-coded. Note that $N$ is a parameter that cannot vary during the continuation. Therefore it does not have entries in Jacobianp. We should let the pde_1 curve know that bruss.m is the active file, the initial values of the parameters $N$ and $L$ are respectively $20$ and $0.06$ and the active parameter is $L$, i.e. the second parameter of bruss.m. So, first of all we have to get the approximate equilibrium solution which is provided in a subfunction 'init' of 'bruss.m'. We first locate the handle to the subfunction init and call it.
>N = 20; L = 0.06;
>handles = feval(@bruss); 
>[t,x0,options] = feval(handles{1},N);
It sets the number of state variables to $2N$ and makes an initial vector $x0$ of length $2N$ containing the values of the approximate equilibrium solution. Now we inform the pde_1 curve that the second parameter of bruss.m is the active parameter and what the default values of the other parameters are. We also set some options.
>[x1,v1] = init_EP_EP(@bruss,x0,[N;L], [2]);
>opt = contset;opt=contset(opt,'MinStepsize', 1e-5);
>opt=contset(opt,'MaxCorrIters', 10);
>opt=contset(opt,'MaxNewtonIters', 20);
>opt=contset(opt,'FunTolerance', 1e-3);
>opt=contset(opt,'Singularities',1);
>opt=contset(opt,'MaxNumPoints',500);
>opt=contset(opt,'Locators',[]);
We start the continuation process by the statement
[x,v,s,h] = cont(@pde_1,x1,v1,opt).
In this case the number of state variables can be a parameter and the Jacobian can be sparse. Using the command
cpl(x,v,s,[41;20]);
we can plot the 20th component of $x$ as a function of $L$ (which is the 41th continuation variable), see Figure 35.

Figure 35: An quilibrium curve of bruss.m
\includegraphics[scale=0.5]{ex/bruss2.eps}

The above test can be executed by running testbrusselator; this file is in the directory Testruns.