ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL

Unit 24 -- Application: Chemical Reactions

Outline of Unit 24

Initialization

> restart;

> with( DEtools ):

> with( plots ):

> with( linalg ):

Warning, the name changecoords has been redefined

Warning, the name adjoint has been redefined

Warning, the protected names norm and trace have been redefined and unprotected

>

An example of a first-order reaction is the conversion of -butyl chloride into -butyl alcohol:

(CH ) CCl + NaOH -> (CH ) COH + NaCl

A chemical reaction is first-order if the molecules of a substance decompose into smaller molecules at a rate proportional to the amount of substance remaining at any time. That is, if denotes the amount of substance at time , then the amount of substance A can be modeled with the first-order linear (and separable) ODE

> ode1 := diff( A(t), t ) = -k * A(t);

>

where the reaction rate is positive and has units of 1/time.

The solution to this ODE, with initial condition

> ic1 := A(0) = A[0];

>

can be written down on sight, but it is simpler to use Maple to enter the result in a Maple worksheet

> sol1 := dsolve( { ode1, ic1 }, A(t) );

>

While first-order chemical reactions are easy to model and analyze, they are not very common.

>

An example of a second-order chemical reaction is

CH Cl + NaOH -> CH OH + NaCl

in which one molecule of methyl alcohol and one molecule of sodium hydroxide combine to form one molecule of methyl hydroxide and one molecule of sodium chloride. This reaction proceeds at a rate proportional to the product of the remaining concentrations of methyl chloride and sodium hydroxide.

The general form for a second-order reaction is

A + B -> C

To model a second-order reaction, let and denote the initial amounts of chemicals and and denote the amount of chemical at time by . The amount of chemical remaining at time is ; likewise, is the remaining amount of chemical . Hence, the amount of C is modeled by

> ode2 := diff( C(t), t ) = k * ( alpha - C(t) ) * ( beta - C(t) );

>

where is a reaction constant. This first-order ODE is nonlinear (but separable). With initial condition

> ic2 := C(0) = C[0];

>

the amount of chemical is found to be

> sol2 := dsolve( { ode2, ic2 }, C(t) );

>

provided . If , the solution becomes

> simplify( dsolve( { eval( ode2, beta=alpha ), ic2 }, C(t) ) );

>

Note that it is not possible to obtain this solution simply by substituting into the general solution

> eval( sol2, beta=alpha );

Error, division by zero

>

but is obtained in the limit as approaches

> map( limit, sol2, beta=alpha );

>

Two molecules of nitrous oxide combine with one molecule of oxygen to form two molecules of nitrogen dioxide:

2NO + O -> 2NO

At room temperature this reaction can be modeled with the IVP

> ode3 := diff( NO2(t), t ) = k * ( alpha - NO2(t) )^2 * ( beta - NO2(t)/2 );

> ic3 := NO2(0) = A;

>

where is the concentration of nitrogen dioxide at time , is the initial concentration of nitrogen oxide, is the initial concentration of oxygen, and is the initial concentration of nitrogen dioxide.

The solution to this IVP with a separable ODE is

> sol3 := dsolve( { ode3, ic3 }, NO2(t) );

>

This "explicit" solution is of almost no use.

Note, however, that the qualitative behavior of solutions can be determined directly from the ODE. There are two equilibria. The equilibrium at is stable and, because of the quadratic term, is semi-stable.

Assuming numerical values for the parameters are available, a numerical solution is easy to obtain.

At a temperature of 25C, the rate constant is 7130 liter /(mole second). Assume the initial concentration of NO is mole/liter, the initial concentration of O is mole/liter and the initial concentration of NO is mole/liter.

> param3 := { k=7.13*10^3, alpha=0.003, beta=0.002, A=0 };

>

A plot of the numerical solution to this IVP for the first hour (3600 seconds) of the reaction is obtained with

> sol3n := dsolve( eval( { ode3, ic3 }, param3 ), NO2(t), numeric ):

> odeplot( sol3n, [t,NO2(t)], 0..3600, view=[0..3600,0..0.003] );

>

Note that the concentration is increasing to the equilibrium solution .

Slightly more information can be obtained by superimposing the solution curve on the direction field for this problem.

> DEplot( eval(ode3,param3), NO2(t), t=0..3600, [[0,0]], NO2=0..0.01, stepsize=25 );

>

Sometimes the slopes are so small that the direction field is not very useful. In such cases, it can be useful to plot the right-hand side of the (autonomous) ODE as a function of the dependent variable.

> plot( eval(rhs(ode3),param3 union {NO2(t)=NO2}), NO2=0.001..0.005, labels=[NO2,"NO2 ' "] );

>

From this plot it is seen that the rate of change of the concentration of NO is positive when the concentration of NO is in [0,0.003) or (0.003,0.004) and the concentration of NO decreases only when the concentration exceeds 0.004.

>