\( \DeclareMathOperator{\abs}{abs} \newcommand{\ensuremath}[1]{\mbox{$#1$}} \)
(%i1) kill(all)$

part 1 : intro staged-process models
[email protected] , 31/01/2017
rem : see : 71e) Lake example. 'second order difference equation'

ricatti equation :     ?  Staged-process

Continuous counter-current Liquid-Liquid extraction

(%i1) batch(solve_rec)$
\[\mbox{}\\\mbox{read and interpret file: C:\ensuremath{\backslash}maxima-5.38.1\ensuremath{\backslash}share\ensuremath{\backslash}maxima\ensuremath{\backslash}5.38.1\_ 5\_ gdf93b7b\_ dirty\ensuremath{\backslash}share\ensuremath{\backslash}solve\_ rec\ensuremath{\backslash}solve\_ rec.mac}\mbox{}\\\mbox{(\%{}i2) eval\_ when(batch,ttyoff:true,nolabels:true)}\]

Figure 1:
Diagram

select : block 'i' → 1,2,..,i,..,n
1e) condition : accumulation = 0
2e)             generation   = 0
3e) phases : complete immiscible (L=water (kg/s),V=kerosene (kg/s))
   in-out + generation = accumulation
   in-out = 0
⇒ (1) material balance 'stag i'
      L*X(i-1)+V*Y(i+1)-(L*X(i)+V*Y(i))=0
⇒ (2) equilibrium relation
      Y(i) = K*X(i)

(%i72) eq1:L*x[i-1]+V*y[i+1]-L*x[i]-V*y[i]=0;
eq2:x[i]=y[i]/k;
\[\tag{eq1}\label{eq1}V\,{{y}_{i+1}}-V\,{{y}_{i}}-L\,{{x}_{i}}+L\,{{x}_{i-1}}=0\] \[\tag{eq2}\label{eq2}{{x}_{i}}=\frac{{{y}_{i}}}{k}\]
(%i74) eq3:subst(y[i]/k,x[i],eq1);
eq4:subst(y[i-1]/k,x[i-1],eq3);
\[\tag{eq3}\label{eq3}-\frac{L\,{{y}_{i}}}{k}+V\,{{y}_{i+1}}-V\,{{y}_{i}}+L\,{{x}_{i-1}}=0\] \[\tag{eq4}\label{eq4}-\frac{L\,{{y}_{i}}}{k}+\frac{L\,{{y}_{i-1}}}{k}+V\,{{y}_{i+1}}-V\,{{y}_{i}}=0\]

def : theta = L/(V*k)

(%i75) eq5:ratsubst(theta,L/(V*k),eq4);
\[\tag{eq5}\label{eq5}V\,{{y}_{i-1}}\theta+V\,{{y}_{i}}\,\left( -\theta-1\right) +V\,{{y}_{i+1}}=0\]
(%i76) eq6:factor(eq5)*(-1);
\[\tag{eq6}\label{eq6}V\,\left( {{y}_{i}}\theta-{{y}_{i-1}}\theta-{{y}_{i+1}}+{{y}_{i}}\right) =0\]
(%i77) eq7:part(eq6,1,2)*(-1)=0;
\[\tag{eq7}\label{eq7}-{{y}_{i}}\theta+{{y}_{i-1}}\theta+{{y}_{i+1}}-{{y}_{i}}=0\]
(%i80) p1:coeff(part(eq7,1),y[i+1]);
p2:coeff(part(eq7,1),y[i]);
p3:coeff(part(eq7,1),y[i-1]);
\[\tag{p1}\label{p1}1\] \[\tag{p2}\label{p2}-\theta-1\] \[\tag{p3}\label{p3}\theta\]

second order : difference equation

(%i81) eq8:p3*y[i-1]+p2*y[i]+p1*y[i+1]=0;
\[\tag{eq8}\label{eq8}{{y}_{i-1}}\theta+{{y}_{i}}\,\left( -\theta-1\right) +{{y}_{i+1}}=0\]

y             = depent      variable
subscript 'i' = independent variable

(%i83) eq9:subst(i+2,i+1,eq8);
eq10:subst(i+1,i,eq8);
\[\tag{eq9}\label{eq9}{{y}_{i-1}}\theta+{{y}_{i}}\,\left( -\theta-1\right) +{{y}_{i+2}}=0\] \[\tag{eq10}\label{eq10}{{y}_{i}}\theta+{{y}_{i+1}}\,\left( -\theta-1\right) +{{y}_{i+2}}=0\]

http://maxima.sourceforge.net/docs/tutorial/en/gaertner-tutorial-revision/Pages/DiffEq0001.htm

(%i84) eq11:solve_rec (eq10, y[i]);
\[\tag{eq11}\label{eq11}{{y}_{i}}={{\mathit{\%{}k}}_{1}}\,{{\theta}^{i}}+{{\mathit{\%{}k}}_{2}}\]

two values : quadratic equation
1e) theta  ⇒ (theta)^i
2e) 1      ⇒ (1)^i

n = integer number ( example n=10)

(%i85) ev(eq11,i=10);
\[\tag{\%{}o85}\label{o85} {{y}_{10}}={{\mathit{\%{}k}}_{1}}\,{{\theta}^{10}}+{{\mathit{\%{}k}}_{2}}\]

How to find :  %k1,%k2
input  x[0]    = > from equilibrium equation 'see above'
output y[n+1]

(%i86) eq12:y[i]:=k*x[i];
\[\tag{eq12}\label{eq12}{{y}_{i}}:=k\,{{x}_{i}}\]
(%i88) eq13:y[0]=subst(0,i,rhs(eq11));
eq14:subst(n+1,i,lhs(eq11))=subst(n+1,i,rhs(eq11));
\[\tag{eq13}\label{eq13}{{x}_{0}}k={{\mathit{\%{}k}}_{2}}+{{\mathit{\%{}k}}_{1}}\] \[\tag{eq14}\label{eq14}{{y}_{n+1}}={{\mathit{\%{}k}}_{1}}\,{{\theta}^{n+1}}+{{\mathit{\%{}k}}_{2}}\]
(%i89) eq15:solve([eq13,eq14],[%k[1],%k[2]]);
\[\tag{eq15}\label{eq15}[[{{\mathit{\%{}k}}_{1}}=-\frac{k\,\left( {{x}_{0}}-{{x}_{n+1}}\right) }{{{\theta}^{n+1}}-1},{{\mathit{\%{}k}}_{2}}=\frac{k\,\left( {{x}_{0}}\,{{\theta}^{n+1}}-{{x}_{n+1}}\right) }{{{\theta}^{n+1}}-1}]]\]
(%i91) eq16:part(eq15,1,1);
eq17:part(eq15,1,2);
\[\tag{eq16}\label{eq16}{{\mathit{\%{}k}}_{1}}=-\frac{k\,\left( {{x}_{0}}-{{x}_{n+1}}\right) }{{{\theta}^{n+1}}-1}\] \[\tag{eq17}\label{eq17}{{\mathit{\%{}k}}_{2}}=\frac{k\,\left( {{x}_{0}}\,{{\theta}^{n+1}}-{{x}_{n+1}}\right) }{{{\theta}^{n+1}}-1}\]
(%i92) eq18:ev(eq11,eq16,eq17);
\[\tag{eq18}\label{eq18}{{x}_{i}}k=\frac{k\,\left( {{x}_{0}}\,{{\theta}^{n+1}}-{{x}_{n+1}}\right) }{{{\theta}^{n+1}}-1}-\frac{k\,\left( {{x}_{0}}-{{x}_{n+1}}\right) \,{{\theta}^{i}}}{{{\theta}^{n+1}}-1}\]
(%i93) eq19:rhs(eq18);
\[\tag{eq19}\label{eq19}\frac{k\,\left( {{x}_{0}}\,{{\theta}^{n+1}}-{{x}_{n+1}}\right) }{{{\theta}^{n+1}}-1}-\frac{k\,\left( {{x}_{0}}-{{x}_{n+1}}\right) \,{{\theta}^{i}}}{{{\theta}^{n+1}}-1}\]
(%i94) eq20:y[i]=eq19;
\[\tag{eq20}\label{eq20}{{x}_{i}}k=\frac{k\,\left( {{x}_{0}}\,{{\theta}^{n+1}}-{{x}_{n+1}}\right) }{{{\theta}^{n+1}}-1}-\frac{k\,\left( {{x}_{0}}-{{x}_{n+1}}\right) \,{{\theta}^{i}}}{{{\theta}^{n+1}}-1}\]

? I don't know how kill the definition,extra steps need it now.
y[i]:= k*x[i]

(%i96) eq21:ratsubst(ti,x[i]*k,eq20);
eq22:ratsubst(tn,x[n+1]*k,eq21);
\[\tag{eq21}\label{eq21}\mathit{ti}=\frac{{{x}_{0}}k\,{{\theta}^{n+1}}+\left( k\,{{x}_{n+1}}-{{x}_{0}}k\right) \,{{\theta}^{i}}-k\,{{x}_{n+1}}}{{{\theta}^{n+1}}-1}\] \[\tag{eq22}\label{eq22}\mathit{ti}=\frac{{{\theta}^{i}}\,\left( \mathit{tn}-{{x}_{0}}k\right) -\mathit{tn}+{{x}_{0}}k\,{{\theta}^{n+1}}}{{{\theta}^{n+1}}-1}\]
(%i98) eq23:subst(yy[i],ti,eq22);
eq24:subst(yy[n+1],tn,eq23);
\[\tag{eq23}\label{eq23}{{\mathit{yy}}_{i}}=\frac{{{\theta}^{i}}\,\left( \mathit{tn}-{{x}_{0}}k\right) -\mathit{tn}+{{x}_{0}}k\,{{\theta}^{n+1}}}{{{\theta}^{n+1}}-1}\] \[\tag{eq24}\label{eq24}{{\mathit{yy}}_{i}}=\frac{{{x}_{0}}k\,{{\theta}^{n+1}}+\left( {{\mathit{yy}}_{n+1}}-{{x}_{0}}k\right) \,{{\theta}^{i}}-{{\mathit{yy}}_{n+1}}}{{{\theta}^{n+1}}-1}\]

i=1 :exit composition 'yy[1]=y[1]'

(%i99) eq25:subst(1,i,eq24);
\[\tag{eq25}\label{eq25}{{\mathit{yy}}_{1}}=\frac{{{x}_{0}}k\,{{\theta}^{n+1}}+\left( {{\mathit{yy}}_{n+1}}-{{x}_{0}}k\right) \theta-{{\mathit{yy}}_{n+1}}}{{{\theta}^{n+1}}-1}\]
(%i100) eq26:solve(eq25,(theta)^(n+1));
\[\tag{eq26}\label{eq26}[{{\theta}^{n+1}}=-\frac{\left( {{\mathit{yy}}_{n+1}}-{{x}_{0}}k\right) \theta-{{\mathit{yy}}_{n+1}}+{{\mathit{yy}}_{1}}}{{{x}_{0}}k-{{\mathit{yy}}_{1}}}]\]
(%i101) eq27:part(eq26,1);
\[\tag{eq27}\label{eq27}{{\theta}^{n+1}}=-\frac{\left( {{\mathit{yy}}_{n+1}}-{{x}_{0}}k\right) \theta-{{\mathit{yy}}_{n+1}}+{{\mathit{yy}}_{1}}}{{{x}_{0}}k-{{\mathit{yy}}_{1}}}\]

how to find :n+1 ,use log on both side 'number of stages'
rem  :    ln(10)=log(10.0)        ' Maxima'
rem  :  def:  log a^n = n*log(a)

(%i102) log(10.0);
\[\tag{\%{}o102}\label{o102} 2.302585092994046\]
(%i103) eq28:phi=rhs(eq27);
\[\tag{eq28}\label{eq28}\phi=-\frac{\left( {{\mathit{yy}}_{n+1}}-{{x}_{0}}k\right) \theta-{{\mathit{yy}}_{n+1}}+{{\mathit{yy}}_{1}}}{{{x}_{0}}k-{{\mathit{yy}}_{1}}}\]

rem : y[i]=yy[i] ,i=1,n+1

(%i104) eq29: n+1 = log(phi)/log(theta);
\[\tag{eq29}\label{eq29}n+1=\frac{\log{\left( \phi\right) }}{\log{\left( \theta\right) }}\]

? How to find operating line : yy(n+1) = a*xx(n)+b
 find :a( = slope),b
see : fig 'above',use for Graphical stage calculations (x[i],y[i]).
     x[i],y[i] : 0 = < x[i],y[i] = < 1
1e) line : operating line
2e) line : equilibrium line
find : number of stages .

in  arrow = (+) sign
out arrow = (-) sign

material balance : stage 1 = L*xx[0]-V*yy[1]
                  stage n = -L*xx[n]+V*yy[n+1]
                  stage 1 = stage n

(%i107) op1:stage1=L*xx[0]-V*yy[1]$
op2:stagen=-L*xx[n]+V*yy[n+1]$
op: rhs(op2)=rhs(op1);
\[\tag{op}\label{op}V\,{{\mathit{yy}}_{n+1}}-L\,{{\mathit{xx}}_{n}}={{\mathit{xx}}_{0}}L-{{\mathit{yy}}_{1}}V\]
(%i108) op3:solve(op,yy[n+1])[1];
\[\tag{op3}\label{op3}{{\mathit{yy}}_{n+1}}=\frac{L\,{{\mathit{xx}}_{n}}-{{\mathit{yy}}_{1}}V+{{\mathit{xx}}_{0}}L}{V}\]

rem yy[n+1]= op4/op5

(%i112) op4:part(rhs(op3),1)$
op5:part(rhs(op3),2)$
op6: a=coeff(op4,xx[n])/op5;
op7: b=ev(op4,xx[n]=0)/op5;
\[\tag{op6}\label{op6}a=\frac{L}{V}\] \[\tag{op7}\label{op7}b=\frac{{{\mathit{xx}}_{0}}L-{{\mathit{yy}}_{1}}V}{V}\]

operating line :

(%i113) op8:yy[n+1]=rhs(op6)*xx[n]+rhs(op7);
\[\tag{op8}\label{op8}{{\mathit{yy}}_{n+1}}=\frac{L\,{{\mathit{xx}}_{n}}}{V}+\frac{{{\mathit{xx}}_{0}}L-{{\mathit{yy}}_{1}}V}{V}\]
Created with wxMaxima.