## Nonfirstorder systems

When a reaction scheme contains second- and higher-order reactions (see Section 1.3.4) it turns out that except in very special cases, there is no general analytical solution to the set of simultaneous differential equations that describe the system. When a differential equation contains the product of two or more dependent variables (concentrations of reactants in the present context) it is said to be nonlinear. For example, the next, relatively simple, reaction scheme cannot be solved analytically because two of its differential equations are nonlinear:

k1 k2

However, the equations can be solved by using a procedure called numerical integration. This is the topic of the next section. In order to provide some motivation to proceed with the present topic it is encouraging to note that in Mathematica, numerical integration of an array of linear or nonlinear differential equations is done by simply adding the letter "N" to the DSolve function, viz., NDSolve.

NDSolve[eqns , find numerical solutions for several functions

{y1, y2, ...}, {x, xmin, xmax}] yi with x in the range xmin to xmax

Finding numerical solutions to differential equations.

Q: Simulate the time course of the reaction scheme in Eqn [1.20] using NDSolve.

A: When performing numerical integration, the output is given as numerical values, so all initial conditions and parameter values must be specified before executing the function or algorithm. There can be no parameter with an unassigned value and this includes specifying the minimum and maximum values of the independent variable, time.

solution = NDSolve[{

a' [t] &-k1a[t] b[t] + k-1c[t], b'[t] &-k1a[t] b[t] + k-1 c[t] , c'[t] & k1a[t] b[t] - k-1c[t] - k2c[t], d'[t] & k2 c[t] , a & a0, b & b0, c & 0, d & 0}, {a[t] , b[t] , c[t] , d[t]} , {t, 0, 10}]

{{a[t] ^InterpolatingFunction[{{0., 10.}}, <>][t], b[t] ^InterpolatingFunction[{{0., 10.}}, <>][t], c[t] ^InterpolatingFunction[{{0., 10.}}, <>][t], d[t] ^InterpolatingFunction[{{0., 10.}}, <>][t]}}

Notice that unlike DSolve, where the solutions for each reactant were given as analytical functions, the solutions for the reactants resulting from NDSolve are given as objects called Interpolating Functions. The Interpolating

Function effectively stores a table of values for each reactant as a function of t, then interpolates in this table to find an approximation to each reactant at the particular t that is requested.

We can plot the results in exactly the same way as we did for DSolve.

Plot[Evaluate[{a[t] ,b[t] , c[t], d[t]} /. solution], {t, 0, 10}, AxesLabel -> {"Time ", "Concentration"}];

Concentration Time

Figure 1.7. Simulation of the time course of the reaction described by Eqn [1.20], using NDSolve operating on the corresponding differential rate equations. The ordinate denotes concentration in mmol L-1 and the abscissa, time in seconds.

Time

Figure 1.7. Simulation of the time course of the reaction described by Eqn [1.20], using NDSolve operating on the corresponding differential rate equations. The ordinate denotes concentration in mmol L-1 and the abscissa, time in seconds.