kill(all)$

part 1 : intro staged-process models
peter.vlasschaert@gmail.com , 31/01/2017
rem : see : 71e) Lake example. 'second order difference equation'

ricatti equation :     ?  Staged-process

Continuous counter-current Liquid-Liquid extraction 

batch(solve_rec)$

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)

eq1:L*x[i-1]+V*y[i+1]-L*x[i]-V*y[i]=0;
eq2:x[i]=y[i]/k;

eq3:subst(y[i]/k,x[i],eq1);
eq4:subst(y[i-1]/k,x[i-1],eq3);

def : theta = L/(V*k)

eq5:ratsubst(theta,L/(V*k),eq4);

eq6:factor(eq5)*(-1);

eq7:part(eq6,1,2)*(-1)=0;

p1:coeff(part(eq7,1),y[i+1]);
p2:coeff(part(eq7,1),y[i]);
p3:coeff(part(eq7,1),y[i-1]);

second order : difference equation 

eq8:p3*y[i-1]+p2*y[i]+p1*y[i+1]=0;

y             = depent      variable 
subscript 'i' = independent variable

eq9:subst(i+2,i+1,eq8);
eq10:subst(i+1,i,eq8);

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

eq11:solve_rec (eq10, y[i]);

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

n = integer number ( example n=10)

ev(eq11,i=10);

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

eq12:y[i]:=k*x[i];

eq13:y[0]=subst(0,i,rhs(eq11));
eq14:subst(n+1,i,lhs(eq11))=subst(n+1,i,rhs(eq11));

eq15:solve([eq13,eq14],[%k[1],%k[2]]);

eq16:part(eq15,1,1);
eq17:part(eq15,1,2);

eq18:ev(eq11,eq16,eq17);

eq19:rhs(eq18);

eq20:y[i]=eq19;

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

eq21:ratsubst(ti,x[i]*k,eq20);
eq22:ratsubst(tn,x[n+1]*k,eq21);

eq23:subst(yy[i],ti,eq22);
eq24:subst(yy[n+1],tn,eq23);

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

eq25:subst(1,i,eq24);

eq26:solve(eq25,(theta)^(n+1));

eq27:part(eq26,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)      

log(10.0);

eq28:phi=rhs(eq27);

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

eq29: n+1 = log(phi)/log(theta);

? 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                

op1:stage1=L*xx[0]-V*yy[1]$
op2:stagen=-L*xx[n]+V*yy[n+1]$
op: rhs(op2)=rhs(op1);

op3:solve(op,yy[n+1])[1];

rem yy[n+1]= op4/op5

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;

operating line :

op8:yy[n+1]=rhs(op6)*xx[n]+rhs(op7);