Solve below in Wolfram Mathematica

a. Prigogine and Lefever described a reaction process that led to oscillations. It is informally referred to as the "Brusselator". The reaction mechanisms are the following: (1) A Overscript[->, Subscript[k, 1]] X;

(2) B + X Overscript[->, Subscript[k, 2]] Y + D;

(3) 2 X + Y Overscript[->, Subscript[k, 3]]3 X;

(4) X Overscript[->, Subscript[k, 4]] E.

Using the law of mass action (see Chapter 5, section 5.6, for a worked example), derive the differential equations for the system. Assume that A and B are in great excess, and can thus be treated as constants; also note that D and E do not generate any products, so that the only differential equations that you need are for the rate of change of X and Y with time.

b. It is possible to nondimensionalize the equations from part a into the following simpler form:

du/d\[Tau]= 1 - (b +1) u + a u^2v (1)

dv/d\[Tau]= b u - a u^2v (2)

where \[Tau] = Subscript[k, 4]t, a = (Subscript[k, 3] (Subscript[k, 1]A)^2)/Subscript[k, 4]^3, and b = (Subscript[k, 2]B)/Subscript[k, 4], u corresponds to X and v corresponds to Y. You may find it helpful to compare the nondimensionalized equations with those you derived in part a. If they differ very significantly, you may want to re-think your derivation in part a. Those of you who are interested are encouraged to show algebraically that these substitutions, and the chain rule, ensure that equations (1) and (2) are dimensionless. It is not necessary to do this, however, to complete the question. In one cell, define a variable parameters, to which you should assign a list of rules for defining the parameters a and b (an example of how to do this is shown in section 5.7 of Chapter 5). Set a to 0.5, and set b to 1.5. In a second cell, define a variable solution, to which you should assign the solution of the two differential equations (1) and (2) evaluated using NDSolve, as well as two additional equations for the initial values for state variables u and v at time 0, which you should set to 0.0. Also make sure to replace the symbolic values of parameters a and b in the equations with those in the variable parameters that you just created using the ReplaceAll operator. Run the system for 200 time units. An example of how to set up a system of two differential equations in NDSolve is also shown in section 5.7 of Chapter 5. In a third cell, use Plot and Evaluate to plot the outputs of the u and v variables over the 200 time units. An example of the use of these commands is shown in section 5.7 of Chapter 5. Describe the results you observe.

c. To explore a much wider range of initial conditions and parameter values, use the commands from part b of this question, suitably modified, as the basis for creating a dynamic model in the function Manipulate (it may help you to review how to do this by looking at Chapter 5, question 2). Set up sliders for the two initial conditions: define variables uinit and vinit for this purpose; set their default values to 1.0 and 3.0, respectively, and give them the titles "Initial value of u" and "Initial values of v", respectively, and ranges 0 to 3. Set up sliders for the two parameters: define variables avalue and bvalue for this purpose; set their default values to 0.5 and 1.5, respectively, their names to "a" and "b", and their ranges from 0.1 to 2 (for the a parameter) and 0 to 2 (for the b parameter). In the plot command, set the plot range to be All for the x axis, and 0 to 25 for the y axis. Also set the PlotPoints option to 200 to ensure that no points are left out as the fastest changing parts of the function are plotted. What happens to the possibility of oscillations, and the amplitude of the oscillations as you vary the initial conditions, or vary a and b? If you need to set any of the parameters to a specific value, remember that you can click on the + sign to the right of the slider for that parameter, and a new line will display, including a box into which you can type exact values.

d. In addition to running simulations, you have a range of powerful analytical tools at your disposal for understanding this system. The first step is to find the equilibrium points for the system. Do this by setting the right hand sides of the differential equations (1) and (2) to zero and solving simultaneously for u and v. Please do this symbolically; you can substitute numbers in for answering the rest of this part, but you will need the symbolic solution below. Explain the result you've obtained. In the model you created in Manipulate in part c of this problem, set parameter a to 0.1, and parameter b to 0.2. Using your result, change the initial condition values so that the outputs u and v do not change over time (i.e., set the system up at an equilibrium point). Describe the results you obtain.

e. To explore the global dynamics of the model qualitatively, we can take advantage of the tools that you have already used: the vector field, nullclines and phase plane plots. Vector field: Define functions ucoord and vcoord whose arguments are u, v, a and b, which are equal to the right hand sides of the differential equations for u and v, respectively. Define a graphics object, gvectorfield, which you should set equal to the vector field plot (using VectorPlot) of ucoord and vcoord, setting arguments a and b to 0.2 and 1.3, respectively, and with u and v both ranging from 0 to 15. What does the flow look like for these settings? Nullclines: To obtain the nullcline for the u variable, set the right hand side of the differential equation to zero, and then use Solve to solve the equation for variable v. Define a function, nullclineu, whose arguments are u, a and b, and set it equal to the solution for the u nullcline. Obtain the v nullcline by setting the right hand side of the differential equation for v equal to zero, and then using Solve to solve the equation for the variable v. Define a function, nullclinev, whose arguments are u, a and b, and set it equal to the solution for the v nullcline. Now create a new graphics variable, gnullclines, and set it equal to a plot of the two nullcline functions. Set a to 0.2, b to 1.3, and have u range from 0.43 to 15. Set the plot range from 0 to 15 for both the x and y axes. Describe the results. Locate the equilibrium point(s) for the system. Use the Show command to superimpose the vector field and nullclines. What can you say about the relationship of the vector field and the nullclines? Trajectories in the phase plane: To see a specific trajectory in the phase plane, run the system as you did in part b of this problem for 200 time units. Set the value of a to 0.2, the value of b to 1.3, the initial condition of u to 1.0, and the initial condition of v to 5.5. Instead of creating a time plot of the result, however, define a new variable, gtrajectory, and set it equal to a parametric plot of the solution (using ParametricPlot). Set PlotPoints to 200, the plot range to 0 to 15 for both the x and y axes, and PlotStyle to Red, so that the trajectory will be plotted in red. Describe the result that you obtain. Using the Show command, superimpose the vector field, the nullclines, and the trajectory. How does the trajectory relate to the vectors in the vector field and to the nullclines? Putting it all together: Create a dynamic model using the pieces you've created so far in this question. In an initialization cell, put the definitions of the functions ucoord, vcoord, nullclineu and nullclinev. Within Manipulate, put in gvectorfield, gnullclines, parameters, solution, gtrajectory, and the Show command to show the vector field, nullclines and trajectory superimposed. Define four variables for manipulate - uinit, vinit, avalue, and bvalue - and modify the commands within manipulate to replace specific values for the parameters and initial values with these variables. With the manipulate command region, set the initial value of uinit to 1.0, name it "Initial u value", and let it range from 0 to 5. Set the initial value of vinit to 5.5, name it "Initial v value", and let it range from 0 to 6. Set the initial value of avalue to 0.2, name it "a", and let it range from 0.1 to 3. Finally, set the initial value of bvalue to 1.3, and let it range from 0 to 2. Play with the a and b parameter values. What happens when b is much larger than a? What happens when b is smaller than a? What qualitative change do you observe under these two different conditions? How do the initial conditions uinit and vinit alter the trajectories under these two different conditions?