% OC.m
% Solves a transient diffusion problem using orthogonal collocation.
% need planar.m

global xoc woc aoc boc qinv nx

nx = 7
planar(nx,1)
for i=2:nx-1
c0(i-1) = 0.;
end

tbreak = [0 .001 .01 .1 1];
for kkk=1:4
tspan = [tbreak(kkk) tbreak(kkk+1)];
[t cc] = ode45('OC_rhs',tspan,c0);
nn=size(t);
c(kkk,:) = cc(nn(1),:);
c0(:) = c(kkk,:);
end

% add in the first and last point
for kkk=1:4
d(kkk,1) = 1;
for i=2:nx-1
d(kkk,i)=c(kkk,i-1);
end
d(kkk,nx) = 0;
end

plot(xoc,d(4,:),'g-*',xoc,d(3,:),'m-*',xoc,d(2,:),'b-*',xoc,d(1,:),'r-*')
xlabel('x')
ylabel('c')
legend('t = 1.0','t = 0.1','t = 0.01','t = 0.001')
title('Diffusion Problem Solved using Orthogonal Collocation')
pause

x = 0:.02:1
for kkk=1:4
pp = spline(xoc,d(kkk,:));
e(kkk,:) = ppval(pp,x);
end

plot(xoc,d(4,:),'g*',xoc,d(3,:),'m*',xoc,d(2,:),'b*',xoc,d(1,:),'r*')
xlabel('x')
ylabel('c')
legend('t = 1.0','t = 0.1','t = 0.01','t = 0.001')
title('Diffusion Problem Solved using Spline Interpolation of Orthogonal Collocation Solution')
hold on
plot(x,e(4,:),'g-',x,e(3,:),'m-',x,e(2,:),'b-',x,e(1,:),'r-')

**********************************************************************
% OC_rhs.m
function ydot=OC_rhs(time,y)
global xoc woc aoc boc qinv nx

% add in the boundary conditions
c(1) = 1;
c(nx) = 0;

% transfer the y to the c
for i=1:nx-2
c(i+1) = y(i);
end

% calculate the right-hand side
for i=2:nx-1
sum = 0;
for j = 1:nx
sum = sum + boc(i,j)*c(j);
end
ydot(i-1) = sum;
end
ydot = ydot';

**********************************************************************
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
if kon==0
fprintf('The total number of collocation points is %g\n\n',nx)
end

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