Problem 26              Due  4/9/01

A model of catalytic converter ignition consists of
the zeroth order kinetic model:

dT/dt = A * exp(-E/kT)

where E is the activation energy, k is Boltzman's
constant, and A is the preexponential rate
constant divided by the heat capacity.  Write an
adaptive algorithm using a local quadrature module
consisting of the G7-K15 quadrature rules to determine
the time necessary to go from a manifold temperature of
600K to the reaction temperature of 900K.  Use a 
tolerance of 1e-8.  Confirm your result by
comparing your algorithm to the matlab routine
'quad'.

Use the constants:
E/k = 9.1x10^3  K
A = 8.9x10^6 K/s

Hint:  obtain the time by integrating the inverse of
the reaction rate.

You will find the following weights and nodes for the
G7-K15 rules useful (based on the domain [-1,1]):

xs=[.949107912342758;.741531185599384;.405845151377397];
xs=[-xs;0;xs];
wxs=[.129484966168870,.279705391489277,.381830050505119];
wxs=[wxs,.417959183673469,wxs];
ys=[.991455371120813;.864864423359769;.586087235467691;...
    .207784955007898];
ys=[-ys;ys];
wxks=[.063092092629979,.140653259715525,.190350578064785];
wxks=[wxks,.209482141084728,wxks];
wyks=[.022935322010529,.104790010322250,.169004726639267,...
      .204432940075298];
wyks=[wyks,wyks];

How to write the quadrature algorithm:  As discussed in class, you
need two parts - the local quadrature module and the control
algorithm.  

1) The local quadrature module takes in the name of the
function to be integrated as a string variable argument, together with
the local range of integration.  It maps the quadrature weights and nodes
to the local domain, and evaluates the named function at the nodes
using the 'feval' command.  It returns the value of the integral over
the domain and an error estimate.

2) The control algorithm takes in the name of the function to be integrated
as a string variable argument, together with the total range of integration
and the tolerance.  Using the 'nargin' function you can actually set the function
up so that it will use some default tolerance if you don't specify one - a convenient
feature.  You need to keep track of the integral, error, and ranges of each sub-
interval as arrays.  If the total error (the sum of the error array) exceeds the
tolerance, you simply find the index corresponding to the maximum error (the
'max' command is useful here!) and subdivide that interval.  I find it simplest
to put the left half in the original spot of the array and to append the right half
to the end of the array.  You then feed these two subintervals to the local
quadrature module and add up the errors again!

Be sure to keep your solution to this problem (your adaptive
quadrature algorithm) in your Matlab toolbox - it's alot
better than the one that Matlab provides!