clear
echo on
%Solution to Problem 26
%
%This problem involved writing an adaptive quadrature
%scheme to evaluate the integral of a function to a
%specified tolerance.  The functions which accomplish
%this are attached below.  The adaptive scheme is given
%in the function myquad.m, the local quadrature module
%is in the file LQM.m, and the function itself (the
%integrand) is in the routine fun.m
pause

%First let's plot up the function:
T=[600:900];
global flag
flag=0; %We add a flag to control a bit of graphics.
figure(1)
plot(T,fun(T))
xlabel('Temperature')
ylabel('dt/dT')
hold on
pause

%Note that most of the time will be taken up going from
%600 to 675 degrees - e.g., where the reaction rate is
%slowest.
pause
%Now we simply call the quadrature function.  We have
%added a bit of graphics to see where the interval is
%being subdivided.
flag=1;
time=myquad('fun',600,900,1e-8);
%The time is thus:
format long
time
hold off
pause

%We compare this to the matlab routine quad.m:
figure(2)
flag=0;
plot(T,fun(T))
xlabel('Temperature')
ylabel('dt/dT')
hold on
flag=1;
quad('fun',600,900,1e-8)

%which is very close.  Note, however, that the quad.m
%rule actually reached its singularity limit (the limit
%of domain discretization), while the G7-K15 rule only
%had to do a single bisection of the domain!  That is
%one reason why the Gaussian quadrature rules are so
%much better!
hold off
echo off



_____________________________


function y=myquad(fname,a,b,tol)

%This function uses the attached local quadrature
%module LQM.m to calculate the integral of the
%function name passed in the string fname from
%a to b with a tolerance of tol.  It uses an adaptive
%algorithm in which the subinterval with the largest
%local error is halved at each step until the error
%tolerance criteria is met.  The function must take
%in an array of nodes and return an array of function
%evaluations as a row vector. If you do not specify
%a tolerance, the function uses a default value of
%1e-4.

if nargin < 4, tol=1.e-4; end
n=1;
alpha(1)=a;
beta(1)=b;
[locint(1),locer(1)]=LQM(fname,a,b);
while sum(locer)>tol*sum(abs(locint)),
 [x,i]=max(locer);
 n=n+1;
 beta(n)=beta(i);
 beta(i)=(alpha(i)+beta(i))/2;
 alpha(n)=beta(i);
 [locint(i),locer(i)]=LQM(fname,alpha(i),beta(i));
 [locint(n),locer(n)]=LQM(fname,alpha(n),beta(n));
end
%And now we have the integral:
y=sum(locint);

___________________________


function [locint,locer]=LQM(fname,alpha,beta)

%This program is a local quadrature module which
%returns the integral of some function given by
%the string 'fname' over the interval [alpha,beta]
%and an estimate of the error.  We employ the G7-K15
%Gaussian quadrature error estimation algorithm.
%The program is called by the command
%[locint,locer]=LQM('f',alpha,beta).

%We set up the Gaussian nodes and weights.  We define
%them on the domain [-1,1].  Note that the order doesn't
%matter as long as we are consistent between weights
%and nodes.
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];

%Now we map to the [alpha,beta] interval:
x=(alpha+beta)/2+(beta-alpha)/2*xs;
y=(alpha+beta)/2+(beta-alpha)/2*ys;
wg=wxs*(beta-alpha)/2;
wk=[wxks,wyks]*(beta-alpha)/2;

%Next we evaluate the function at the node locations:
fx=feval(fname,x);
fy=feval(fname,y);

%And we calculate the quadratures:
g7=wg*fx;
k15=wk*[fx;fy];
%Where:
locint=k15;

%Now for the error estimate.  From Kahaner, Moler
%and Nash we have an error estimate based on the
%difference between the g7 and k15 rules.  We will
%scale this with the magnitude of f:
scale=abs((beta-alpha)*norm(fx));
locer=scale*(200*abs((g7-k15)/scale))^1.5;
%And that's it!



________________________


function y=fun(T)

%This function yields the integrand of problem
%22.  It is the inverse of the rate of change
%of the temperature.  Thus:
global flag

A=8.9e6;
Ek=9.1e3;
y=exp(Ek./T)/A;

if flag==1,plot(T,y,'o');end