kill(all)$

maxima :concentration respons for distilation column ,lsoda 'version 1.0'
peter.vlasschaert@gmail.com,30/01/2021

method : lsoda

see:**
===
------------------------------------------------------------------------------------------------------------------------------------------
maxima :concentration respons for distilation column ,version 1.0
peter.vlasschaert@gmail.com,26/08/2020
------------------------------------------------------------------------------------------------------------------------------------------
------------------------------------------------------------------------------------------------------------------------------------------
https://common-lisp.net/~rtoy/maxima/docref/maxima_75.html#Item_003a-Introduction-to-ODEPACK
------------------------------------------------------------------------------------------------------------------------------------------

load packages & option
===================

load(lrats)$
ratprint :false $

y[i],relative volatility =α
=================
y[i]=f(α,x[i])

m(α,i):=(α*x[i])/(1+(α-1)*x[i]);

y[i], i = 1,2,..4

makelist(y[i]:m(α,i),i,1,3);

y[4]:x[4]; /*total condensor */

differential equations : distillation column (four stage :distillation : column)
========================================================
**

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     */

m1:makelist(rhs(q1[i]),i,1,4);

m2:listofvars(m1);

input : differential for distillation
========================

→ constant

α:4.78 $     /*constant : relative volatility */
F:1$         /* Feed Rate                     */
zf:0.5 $     /*Feed Composition               */

=>  Vapour flow

V:3.55 $ V1:V $ V2:V $ V3:V $

=> Liquid flow

L:3.05 $ L4:L $ L3:L $ L2:L+F $

→  equation  Bottom,Distillate

D : V3-L4$
B : L2-V1$

differential equations with (input)
=========================

m3:ev(m1,nouns);

m4:makelist(x[i]:concat(x,i),i,1,length(m3));

m5:ev(m3,nouns);

m6:listofvars(m5);

start condition : distillation ( ≠ IC intial condition)
====================================

m7:[0.5d0,0.5d0,0.5d0,0.5d0];

1e)  step : load package
----------------------------------

load("dlsode")$

2e)  step :  variables
-----------------------------

m8:flatten([t,m6])$

3e) step :state
---------------------

state : dlsode_init(m5, m8, 21)$

rem: general select : 21, see : state
-----------------------------------------------------------------------
10:Nonstiff (Adams) method, no Jacobian used.     
21:Stiff (BDF) method, user-supplied full Jacobian.     
22:Stiff method, internally generated full Jacobian. 
------------------------------------------------------------------------

4e) step:init values method
--------------------------------------

x :m7$
t: 0d0$
rtol : 1d-6$
atol: [1d-6, 1d-6, 1d-6,1d-6]$ /* length same : length(m4) */
istate: 1$   /* see ref : above */
t:0d0$
stepv:0.1d0$ /* stepvalue */
tout:stepv$  /* d0, double float */

5e) step : simulation method.
----------------------------------------

nn:11$/*number of steps-1 */

m10:flatten([t,m7]);

m:[m10]$ /*init data : from  simulation  */

6e)  step : output simulation.
----------------------------------------

→  init values

printf(true, "At time = ~20,3,2e  x = ~{~14,3,2e~}", 0, m7)$

→  interactive output

for k : 2  thru nn do
  block([],
    result: dlsode_step(x, t, tout, rtol, atol, istate, state),
      q:flatten([result[1],result[2]]),
      m:endcons(q,m),
    printf(true, "At time = ~20,3,2e  x = ~{~14,3,2e~}~%", result[1], result[2]),
    tout : k*stepv);

example : block inside for loop construction. (see:ref)

rem :printf : 20 & 14 you can change without problem for the simulation.

→  fixed output (list,matrix)

list : output
---------------

m11:transpose(m);

matrix :output
--------------------

m12:apply(matrix,m);

how get data matrix: m12
====================

matsize(A):=[length(A),length(transpose(A))];

m12[1];
m12[2,2];

matsize(m12);

plot 2d : discrete, kk=2,3,...length(m4)+1
===============================

→ create data : discrete plot2d

txdata(kk):=makelist([m12[ll][1],m12[ll][kk]],ll,1,matsize(m12)[1]);

→ plot2d

wxplot2d([[discrete,txdata(2)],[discrete,txdata(3)],[discrete,txdata(3)],[discrete,txdata(4)],[discrete,txdata(5)]]);

--------------------------------------------------------------------------------------------------------------------------------------------------------------

first time : use printf '2021', formated interactive output.
---------------------------------------------------------------------------
used iteration output : ~{~14,3,2e~}, x1,x2,x3,x4
---------------------------------------------------------------------------
~% : go to new line
----------------------------------------------------------------------------