clear
echo on
%Solution to problem 5
%This script determines the minimum value
%of chi for which 99% of the inlet concentration
%of solute is removed from the water stream.
%We use the attached function xout.m which
%calculates the outlet concentration
%as a function of chi.  We determine the minimum
%value of chi both graphically and using Matlab's
%root finding command fzero.  Type in 'help fzero'
%to learn more about it.
%Thus:
chi=[0:.1:4];
semilogy(chi,xout(chi))
xlabel('chi')
ylabel('outlet concentration')
grid
pause
%From this plot it is seen that the concentration
%reaches 1% of the inlet concentration at about
%chi=2.8.  We can do better than this using the
%fzero command:
chi_min=fzero('xout(x)-0.01',2.8)
%where 2.8 is our initial guess.  The actual
%value is close to that obtained from looking at the
%graph.
hold on
semilogy(chi_min,xout(chi_min),'o')
hold off
echo off


_________________________


function y=xout(chi)
%This function illustrates how a system of equations
%can be used to solve a chemical engineering problem.
%Consider a liquid-liquid extraction unit.  In each stage
%of the unit the water and the solvent streams are
%assumed to be in equilibrium, such that y(i)=k*x(i)
%where k is the henry's law constant.  In addition,
%we have a mass balance on each stage which basically
%states that the amount of solute going into a stage
%equals to the amount going out - usually a good idea.
%
%We want to use this information to determine the
%concentration distribution in a four stage unit.
%
%We let s be the molar flow rate of the solvent stream
%and we let w be the molar flow rate of the water
%stream.  k is the henry's law constant.  The
%solution depends solely on the ratio s*k/w.  We
%give this parameter the symbol chi.  This is the input
%variable for the function.
%
%We shall adapt the function to handle an array of input
%values of chi.
npt=length(chi);
for kk=1:npt
chit=chi(kk);

%
%We set up the matrix:

n=4;
%This line of code creates the main diagonal of the
%matrix we wish to solve:
a=-(1+chit)*diag(ones(1,n),0);
%
%This line creates the sub-diagonal of the matrix:
a=a+diag(ones(1,n-1),-1);
%
%And this line creates the super diagonal:
a=a+chit*diag(ones(1,n-1),1);
%
%We set up the right hand side:
b=zeros(n,1);
b(1)=-1;
%where we have assumed that the inlet concentration
%of the solute in the water is unity, and that the
%inlet concentration of the solute in the solvent
%is zero.
%
%We now solve for the concentration distribution
%using the \ command:
x=a\b;
%
%and finally, we return the concentration of the nth
%stage:
y(kk)=x(n);
%
%and then we loop back to the next value of chi:
end;