kill(all)$

maxima: matrix structure 'concentration respons dist.column',version 1.0
peter.vlasschaert@gmail.com,26/08/2020

load packages ,option
=================

load(lrats)$ ratprint:false$

see :website
maxima :concentration respons for distilation column ,version 1.0
peter.vlasschaert@gmail.com,26/08/2020

differential equations : distillation column (four stage :distillation : column)
========================================================
unknown : x[i], i=1..nd, matrix structure ( time depent => dynamics)

nd:4 $ /*(n)umber of (d)ifferential equations */;

system (d)
---------------

q1[1]:dxdt(1) = L2*x[2]-V1*y[1]-B*x[1];         /*bottom of column */
q1[2]:dxdt(2) = L3*x[3]+V1*y[1]-L2*x[2]-V2*y[2]+F*zf;  /*feed tray*/
q1[3]:dxdt(3) = L4*x[4]+V2*y[2]-L3*x[3]-V3*y[3];/*above feed tray*/
q1[4]:dxdt(4) = V3*y[3]-L4*x[4]-D*y[3];        /*top of column     */

find : lhs 'system (d)'
----------------------------

p1a:makelist(lhs(q1[i]),i,1,nd );

find : rhs 'system (d)'
----------------------------

p1b:makelist(rhs(q1[i]),i,1,nd );

a) insert equilibrium relation : simple  ' y[i] = K[i]*x[i]
-----------------------------------------------------------------------

p2a:makelist(y[i]=K[i]*x[i],i,1,nd);
p2b:lratsubst(p2a,makelist(q1[i],i,1,nd));

new differential equations  ( =system new d)
-----------------------------------------------------------

qq1[1]:p2b[1];
qq1[2]:p2b[2];
qq1[3]:p2b[3];
qq1[4]:p2b[4];

p2c:sort(listofvars(p2b));

p2c1:makelist(x[i],i,1,nd);

matrix :structure ?
-------------------------

tridiagonale system 'm1'

m1:coefmatrix(p2b,p2c1);

missing component : feed tray info ?  term  zf (known : variable)
rem : z = some value

p2c2:append(p2c1,[p2c[length(p2c)]]);

m2:coefmatrix(p2b,p2c2);

rem :augmented matrix for the steady system 
===================================
only difference : steady system ( = time indepent)
-F => F  (because : rhs of non differential form 'dxdt(i) =0')

to find true augmented matrix : with condition from above
( = time independent system ' = algebraic equations')

1e) step : dxdt(i) =0, i =1,2,..nd, ' new system d'
-----------------------------------------------------------------

p3a:makelist(dxdt(i)=0,i,1,nd);

p3a1:lratsubst(p3a,makelist(qq1[i],i,1,nd));

2e) step :find number algebraic equation for feed tray
-------------------------------------------------------------------------

a)  differentation with respect  zf

p3a2:diff(lhs(p3a1),zf);

b) clean up ( = remove the 'lhs')

p3a3:makelist(rhs(p3a2[i]),i,1,nd);

c) find number of algebraic equation, zf

p3a4:ftray(h):=block(
     h1 : length(h),
    for i:1 thru nd do(
    if h[i]≠ 0 then return(i)
    )
)$

optimized code :

grind(ftray);

p3a5:ftray(p3a3);

3e) modified : algebraic equations of feed tray.
----------------------------------------------------------------
use :'p3a1'

a) find equation ' every algebraic eq. : lhs = 0'

p3a6:p3a1[ftray(p3a3)];

b) final modified equation  
use :'p3a3' ( oposit sign=* (-1) ),'p3a5'

p3a7:p3a3[p3a5]*(-1)*zf;

p3a8:p3a6+p3a7;

4e)  get new list with modified equation.
-------------------------------------------------------
use :'p3a1'

p3a9:makelist(if i = p3a5 then 
    p3a8 else p3a1[i],i,1,nd);

reverse :system : rhs =lhs & lhs = rhs

p3a10:makelist(rhs(p3a9[i])= lhs(p3a9[i]),i,1,nd);

reason : unknown = lhs  , known = rhs( = input feed tray ,values)

5e) final result augmented matrix
----------------------------------------------
use :'p2c1' (a), 'p2c2' (b)

a)  [x[1],x[2],x[3],x[4]]

m3a:augcoefmatrix(p3a10,p2c1);

b)  [x[1],x[2],x[3],x[4],zf]

m3b:augcoefmatrix(p3a10,p2c2);

packages:eigen

load(eigen)$

Function: submatrix
==================================================
    submatrix (i_1, …, i_m, M, j_1, …, j_n)
    submatrix (i_1, …, i_m, M)
    submatrix (M, j_1, …, j_n)
==================================================
    Returns a new matrix composed of the matrix 
M with rows i_1, …, i_m deleted, and columns j_1, …, j_n deleted. 
==================================================
use : load(eigen)
==================================================

'true augmented matrix' : time indepent system of algebraic equations.

m3b1:submatrix(m3b,nd+2);

'see rem : 'm2' '

info :additional : multiplication of matrix with list : operator = .
===============================================

m3b2:m3b1.p2c2;

=>  not any physical meaning :m3b2 (easy to multiply with list)

m3b2a:m3b2[1];

m3b3:m3b1.transpose(p2c2);

is(m3b2=m3b3);

result : use of transpose 'redundant', because of the "true" evaluation.