% run_diffusion.m
% This program runs a code to solve for unsteady diffusion.
% It uses orthogonal collocation on finite elements.

global aoc boc qinv xoc woc np ne delx psilong

% np is the total number of collocation points in one finite element
np = 4;
% ne is the number of elements in the x-direction
ne = 8;
% This can be changed to produce a variable grid
for k=1:ne
delx(k) = 1/ne;
end

% number of interior points is np-2 per element
% this is the number of ODEs being solved
nint=(np-2)*ne;

% create the collocation matrices for the fluid
planar(np,0)

% set up the x values got plotting purposes
num=0;
xstart = 0.;
for k=1:ne
for j=1:np-1
num=num+1;
x(num)=xstart+delx(k)*xoc(j);
end
xstart = xstart+delx(k);
end
x(nint+ne+1) = 1.;

% set the initial conditions
for i=1:nint
c0(i) = 0.;
end

tbreak = [0 .001 .01 .1 1];
for kkk=1:4
tspan = [tbreak(kkk) tbreak(kkk+1)];
[t cc] = ode45('rhs_diffusion',tspan,c0);
nn=size(t);
c(kkk,:) = cc(nn(1),:);
c0(:) = c(kkk,:);
smooth(np,ne,delx,c0)
yplot(kkk,:) = psilong;
tbreak(kkk+1)
end

plot(x,yplot(4,:),'r*-',x,yplot(3,:),'b*-',x,yplot(2,:),'m*-',x,yplot(1,:),'g*-')
xlabel('x')
ylabel('c')
title('Orthogonal Collocation on Finite Elements')
legend('t = 1.0','t = 0.1','t = 0.01','t = 0.001')
**********************************************************************
% rhs_diffusion
% This subroutine computes the rhs for the
% diffusion equation using the orthogonal colloction method.

function cdot=rhs_diffusion(t,c)
global aoc boc qinv xoc woc np ne delx ce psilong

% c are the concentrations at the points internal to the elements

% smooth, using the boundary conditions
smooth(np,ne,delx,c)

% calculate the rhs, using the element form
num=0;
for k=1:ne
for j=2:np-1
% compute the right-hand sides for interior collocation points
sum=0.;
num=num+1;
for i=1:np
sum=sum+boc(j,i)*ce(k,i);
end
cdot(num)=sum/(delx(k)*delx(k));
end
end
cdot = cdot';
**********************************************************************
% this subroutine takes the solution at points interior
% to each element, psi, and the two boundary conditions
% and constructs the values at the points between elements
% and the two ends so that the whole function has continuous
% first derivatives throughout and satisfies the boundary
% conditions.
% The boundary conditions at x = 0 and x = 1
% are specified values (1 and 0, respectively).
% np is the number of collocation points in one element
% ne is the number of elements
% delx is the element size, delx(i),i=1,ne
% c is the function at the points interior to each element,
% with (np-2)*ne values.
% Output is psilong(i), at all fluid collocation points

function smooth(np,ne,delx,psi)
global afd bfd cfd dfd
global xoc woc aoc boc qinvoc ce psilong

nt = ne + 1;

% move the c values to the values for each element
for k = 1:ne
for j=2:np-1
ce(k,j) = psi((np-2)*(k-1)+j-1);
end
end

% set up the matrix
for i=2:nt-1
ratio = delx(i-1)/delx(i);
afd(i)=aoc(np,1);
bfd(i)=aoc(np,np)-aoc(1,1)*ratio;
cfd(i)=-aoc(1,np)*ratio;
s1 = 0;
s2 = 0;
for j=2:np-1
s1 = s1 + aoc(np,j)*ce(i-1,j);
s2 = s2 + aoc(1,j)*ce(i,j);
end
dfd(i)=ratio*s2 - s1;
end
% fix the boundary conditions
bfd(1)=1.;
cfd(1)=0.;
dfd(1)=1.;
afd(nt)=0.;
bfd(nt)=1.;
dfd(nt)=0.;

% solve the tridiagonal system
tridiag(nt)

% move into the element form and one long vector
ce(1,1) = dfd(1);
ce(ne,np) = dfd(nt);
for k=2:nt-1
ce(k,1)=dfd(k);
ce(k-1,np) = dfd(k);
end

% move into one long vector
num=0;
for k=1:ne
for j=1:np-1
num=num+1;
psilong(num)=ce(k,j);
end
end
psilong(num+1)=ce(ne,np);
**********************************************************************
function tridiag(n)
global afd bfd cfd dfd

% this subroutine solves a tridiagonal system of
% equations of the type that often occur with the
% finite difference method.
% input
% a(n), b(n), c(n), d(n) in the following equaiton:
% a(i)*y(i-1) + b(i)*y(i) + c(i)*y(i+1) = d(i)
% n - the size of the system being solved.

% output
% afd(n), bfd(n), cfd(n) this is the lu decomposition.
% the original afd, bfd, cfd are destroyed.
% dfd(n) is the solution.

% lower decomposition tridiagonal matrix
for l=2:n
s = afd(l)/bfd(l-1);
bfd(l) = bfd(l)-s*cfd(l-1);
afd(l) = s;
end

% this is the lower decomposition. the routine could
% be separated here if you have several right-hand sides
% with the same matrix.

% forward sweep tridiagonal matrix
for l=2:n
dfd(l) = dfd(l)-afd(l)*dfd(l-1);
end
% back substitution
dfd(n) = dfd(n)/bfd(n);
for l=2:n
k = n-l+1;
dfd(k) = (dfd(k)-cfd(k)*dfd(k+1))/bfd(k);
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