kill(all)$

appendix :control :model:mixing cstr ' with rk(4): 'numerical', step respons',version 1.0a
peter.vlasschaert@gmail.com,06/01/2020

see:
appendix : electrical  , solve one differential equation , numerical & symbolic,version 1.0

load packages
-----------------------

load(diffeq)$

model :mixing model
================
Caf :feed concentration = 1
y = Ca 'concentration a'
x = time 
q=volumetric flow rate = 100 m^3/hr
V=volume of CSTR      = 100 m^3

rem : CSTR = Continuous Stirred-Tank Reactor

d ( Ca)/dt  = q/V*(Caf -Ca)

maxima : function  , 'rk'
=================

n:10;
h:1;
x0:0;
y0:0;

sol: rk(1-y, y, 0, [x, x0, n, h]);
wxplot2d([discrete, sol])$

discrete 

f(x,y,q,V,Caf) := q/V*(Caf-y);

q:100;
V:100;
Caf:1;

f(x,y,q,V,Caf) ;

'diff(y,x) = f(x,y,q,V,Caf) ;

maxima : function  , 'rk'
=================

n:10;
h:1;
x0:0;
y0:0;

 vx: makelist(0, n + 1);   /* reserved array */
 vy: makelist(0, n + 1);  
   x: x0; 
   y: y0;
   vx[1]: x0;
   vy[1]: y0;
   for i from 1 thru n do (
      k1: ev(h * f(x, y,q,V,Caf),float),
      k2: ev(float(h * f(x + h / 2, y + k1 / 2,q,V,Caf))),
      k3: ev(float(h * f(x + h / 2, y + k2 / 2,q,V,Caf))),
      k4: ev(float(h * f(x + h, y + k3,q,V,Caf))),
      y: y + (k1 + 2 * k2 + 2 * k3 + k4) / 6,
      vy[i + 1]: y,
      x:x+h ,         
      vx[i + 1]: x
                                  )$

vx;
vy;

v:[]$
for i from 1 thru n+1 do(
v:endcons([vx[i],vy[i]],v))$
v;
wxplot2d([discrete,v])$

input respons : Caf = 1(fixed),all time (vx[i])
=================================

vv:[]$
for i from 1 thru n+1 do(
vv:endcons([vx[i],Caf],vv))$
vv;

wxplot2d([discrete,vv])$

values  => removed
===============

x;

y;

kill(x,y);

x;
y;

wxplot2d([[discrete,vv],[discrete,v]],[x,-1,10],[y,0,1.2],
    [legend,"Caf","Ca"],[xlabel,"time"])$

[k1,k2,k3,k4];

kill(k1,k2,k3,k4);

[k1,k2,k3,k4];

g(x,y,q,V,Caf) := q/V*(Caf-y);

rk4(g, a1,  a2, h, n) := block(
       x:a1 ,
       y:a2,                 
      k1: ev(h * g(x, y,q,V,Caf),float),
      k2: ev(float(h * g(x + h / 2, y + k1 / 2,q,V,Caf))),
      k3: ev(float(h * g(x + h / 2, y + k2 / 2,q,V,Caf))),
      k4: ev(float(h * g(x + h, y + k3,q,V,Caf))),
          y: y + (k1 + 2 * k2 + 2 * k3 + k4) / 6,
         vyy[j + 1]:y,
         x:x + h,
         vxx[j + 1]: x    
                                                 )$

rem : no index variable allowed rk4  arguments
rk4(g, vxx[j],  vyy[j], h, n) =>  result  'error'

   
   vxx: makelist(0, n + 1);
   vyy: makelist(0, n + 1);
   vxx[1]: x0;
   vyy[1]: y0;
   for j from 1 thru n  do (
      rk4(g, vxx[j],  vyy[j] , h, n) 
                                     )$

vxx;
vyy;

plot discrete values
===============

1e) unit  respons (&t=0, here x)
-------------------------------------------

vv;

2e) values from vxx & vyy
-------------------------------------------

vvvv:[]$
for i from 1 thru n+1 do(
vvvv:endcons([vxx[i],vyy[i]],vvvv))$
vvvv;

3e)  discrete plot
------------------------

kill(x,y);

[x,y];

rem : x,y symbolic  for wxplot2d option : ranges (x,y)

wxplot2d([[discrete,vv],[discrete,vvvv]],[x,-1,10],[y,0,1.2],
    [legend,"Caf","Ca"],[xlabel,"time"])$