kill(all)$

biotech : advection equation, upwind scheme,version 1.0
peter.vlasschaert@gmail.com,24/05/2020

load  packages
=============

load(lrats)$

ratprint:false$

model : advection (upwind scheme : v=1) ,uu=uu(x,t)
=======================================
∂uu/∂t +v*∂uu/∂x = 0 ,   domain      0 ≤ x ≤ 1.0
=======================================
ic  : u(x,0) = exp((-200)*(x-0.25)^2)
bc : periodic condition , 'see below'
=======================================

p1:'diff(uu,t)+v*'diff(uu,x) =0;

info:https://en.wikipedia.org/wiki/Upwind_scheme

given : parameters
==============

xmin:0;       /*begin x:domain */
xmax:1;     /*end  x:domain */
N:100;       /*number of nodes, x values*/
dt:0.009;   /*nstep = tmax/dt) */
t:0;            /* start time */
tmax:0.5; /*max simulation time => number of steps :code */
v:1;           /*wave speed*/

dx:(xmax-xmin)/N;

find : courant number  <= 1
=====================

courant :v*dt/dx;

initial function: u(x,0) = f(x)
-----------------------------------

x:makelist(xmin+ii*dx,ii,0,N);

length(x);

f(xa):=exp(-200*(xa-0.25)^2);

initial matrix 
=========

u:makelist(f(x[kk]),kk,1,length(x))$
un:makelist(f(x[kk]),kk,1,length(x));

length(u);

simulation
========

→  number of time steps

nstep:round(tmax/dt);

→ data:collection

uu:genmatrix(lambda([i,j], 0), nstep, length(u))$

→ program

for n:1 thru nstep do(
       /* periodic boundary conditions */
        un[1]:u[1]-v*(u [N+1]-u[N]),
        un[N+1] :un[1],
     /* scheme upwind scheme : for v>0 */
      for i:2 thru N+1 do(
       un[i] : u[i]-courant*(u[i]-u[i-1])
                                      ),
    for ii:1 thru length(u) do
    (
      /* old = new */
      u[ii] :un[ii] ,
     /* collection of data */
      uu[n,ii]:u[ii]
    )
)$

→ matrix 'info'

n : number of time   steps
i  : number of space steps

uu$

output : multipleplot, uu
====================

→ find row : i , matrix 'uu'

eq[ii]:=makelist(uu[ii,jj],jj,1,length(un));

→ empty list

m2:[]$

→  moving :wave front( = gauss function)

for ii:1 thru nstep step 7 do
(
                               (                                   
        m2:endcons([discrete,x,eq[ii]],m2)
                                )
 )$
 wxplot2d(m2,[xlabel,"x"],[ylabel,"row output uu=eq[ii]"]);

rem : →  different color  = different time   =>   eq[ii],ii=1...nstep, step =7