% filename axialOC.m
% This code solves the equations for a reactor with
% axial dispersion using the orthogonal collocation method.
% need planar.m, planar_interp.m
% If the problem is linear, set linear = 1 to avoid a second iteration.
global xoc woc aoc boc qinv


% set linear to 1 if problem is linear - this avoids one unnecessary iteration
linear = 0
% set the type of reaction rate
iwhich = 2
%set nx, Peclet (Pe), and Damkohler (Da) before calling
peinv = 1/Pe;
planar(nx,1)

% iterate
tol = 1.e-12;
change = 1.;
iter=0;
% initial guess
for i=1:nx
c(i) = 0.0;
end

while change>tol
iter=iter+1;

% Set up the matrices
% set the matrices to zero
for i=1:nx
f(i) = 0.;
for j=1:nx
aa(i,j) = 0.;
end
end
% set the matrices for the differential equation
for i=2:nx-1
for j=1:nx
aa(i,j) = boc(i,j)*peinv - aoc(i,j);
end
ans = rate(c(i),iwhich);
raterxn = ans(1);
drate = ans(2);
aa(i,i) = aa(i,i) - Da*drate;
f(i) = Da*(raterxn - drate*c(i));
end
for j=1:nx
aa(1,j) = - aoc(1,j)*peinv;
aa(nx,j) = aoc(nx,j);
end
ans = rate(c(1),iwhich);
raterxn = ans(1);
drate = ans(2);
aa(1,1) = aa(1,1) + 1;
f(1) = 1;

% Solve one iteration.
c11 = aa\f';

% Calculate the criterion to stop the iterations.
sum=0.;
for i=1:nx
sum = sum + abs(c11(i) - c(i));
c(i) = c11(i);
end
if linear==1,break,end
change = sum
end

% Plot the solution.
% Find the coefficients
dd1 = qinv*c';
xp = 0:0.02:1; %plotting points
cinterp = planar_interp(nx,dd1,xp,0,1);
plot(xoc,c,'*b',xp,cinterp,'-b')
xlabel('x')
ylabel('c')
c

*********************************************
% rate.m
% calculate the rate and derivative for a given concentration

function ans = rate(c,iwhich)

if iwhich==1
% use for linear reaction
ans(1) = c;
ans(2) = 1.;

elseif iwhich == 2
% use for second order reaction
ans(1) = c*c;
ans(2) = 2.*c;
else
% use for non-isothermal reaction
%ggg = 18.; bbb = 0.3;
ggg = 30; bbb = 0.4;
t = 1 + bbb - bbb*c;
eee = exp(ggg*(1 - 1/t));
ans(1) = c*eee;
ans(2) = eee*(1 - bbb*ggg*c/(t*t));
end

*********************************************
function planar(nx,kon)
% this subroutine computes the matrices for collocation without
% symmetry, using collocation points in table 4-3.
% input variable
% nx = number of collocation points, including two end points
%
% output variables, must be declared global
% aoc = matrix for first derivative, eq. 4-103
% boc = matrix for second derivative, eq. 4-103
% qinv = matrix for q inverse, eq. 4-101
% xoc = vector of collocations points, from table 4-3
% woc = vector of weights, eq. 4-106
%
global xoc woc aoc boc qinv
if(nx<3)
nx=3;
fprintf('\n ***************************************** \n')
fprintf('The number of collocation points must be between \n')
fprintf('3 and 7 (inclusive). nx has been changed to %g\n\n',nx)
end
if(nx>7)
nx=7;
fprintf('\n ***************************************** \n')
fprintf('The number of collocation points must be between \n')
fprintf('3 and 7 (inclusive). nx has been changed to %g\n\n',nx)
end
% now 3<=nx<=7
fprintf('The total number of collocation points is %g\n\n',nx)

xoc(1)=0.;
if (nx==3)
xoc(2)=0.5;
elseif (nx==4)
xoc(2)=0.21132486540519;
xoc(3)=1-xoc(2);
elseif (nx==5)
xoc(2)=0.11270166537926;
xoc(3)=0.5;
xoc(4)=1-xoc(2);
elseif (nx==6)
xoc(2)=0.06943184420297;
xoc(3)=0.33000947820757;
xoc(4)=1-xoc(3);
xoc(5)=1-xoc(2);
else
xoc(2)=0.04691007703067;
xoc(3)=0.23076534494716;
xoc(4)=0.5;
xoc(5)=1-xoc(3);
xoc(6)=1-xoc(2);
end
xoc(nx)=1.;
format long
for i=1:nx,
r(i,i) = 0.0;
aoc(i,i) = 0.0;
s(i) = 1.0;
boc(i,i) = 0.0 ;
for j=1:nx,
if (i~=j)
r(i,j) = 1.0/(xoc(i)-xoc(j)) ;
s(i) = s(i)*r(i,j);
end
end
for j=1:nx,
jx = nx-j+1;
if (jx<j), break, end
if (jx==j)
aoc(i,i) = aoc(i,i)+r(i,j);
end
if (jx>j)
aoc(i,i) = aoc(i,i)+r(i,j)+r(i,jx);
end
end
end
for i=1:nx
for j=1:nx
if (i~=j)
aoc(i,j) = s(j)*r(i,j)/s(i);
boc(i,j) = 2.0*aoc(i,j)*(aoc(i,i)-r(i,j));
boc(i,i) = boc(i,i)+r(i,j)*(aoc(i,i)-r(i,j));
end
end
end
for i=1:nx
qinv(1,i) = s(i);
k = 1;
woc(i) = 0.0 ;
for j=1:nx
if (j~=i)
l = k;
k = k+1;
qinv(k,i) = qinv(l,i);
while (l>1)
m = l-1;
qinv(l,i) = qinv(m,i)-xoc(j)*qinv(l,i);
l = m;
end
qinv(1,i) = -xoc(j)*qinv(1,i);
end
end
for j=1:nx
woc(i) = woc(i)+qinv(j,i)/j;
end
end

if kon == 0
fprintf('A matrix (aoc) \n')
disp(aoc)
fprintf('B matrix (boc) \n')
disp(boc)
fprintf('Q matrix (qinv) \n')
disp(qinv)
fprintf('Collocation points (xoc) \n')
disp(xoc)
fprintf('W matrix (woc) \n')
disp(woc)
end

*********************************************
% filename planar_interp.m
function y = planar_interp(nx,dd,x,x0,x1)
% nx = number of collocation points, including end points
% dd = the coefficients in the polynomial (vector)
% x is a vector of the x values to evaluate the polynomial at
% The points start at x0 and go to x1, and the xb below has to
% be scaled to go from 0 to 1.

xb = (x-x0)/(x1-x0);
nn = size(xb);
for j = 1:nn(2)
x = xb(j);
sum = 0;
for i=1:nx
sum = sum*x + dd(nx-i+1);
end
y(j) = sum;
end