kill(all)$

biotech :numerical simple diffusion pde ,version 1.0
peter.vlasschaert@gmail.com,20/05/2020

reference :  'finite difference'
======================

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

youtube : https://www.youtube.com/watch?v=ZSNl5crAvsw
my website : 
see : SPECIAL TOPICS  :  
          PDE  :    numerical

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

load package
===========

load(lrats)$

numerical :definitions
------------------------------

definition
=======

1e)  def : nuelem
------------------------

nuelem(i,nv,xmin,xmax):=
block(x[i]:ev(xmin+(i-1)*(xmax-xmin)/(nv-1)))$

grind(nuelem)$

2e) def : linspace
-------------------------

linspace(xmin,xmax,nv):=block(
ev(makelist(nuelem(i,nv,xmin,xmax),i,1,nv),nouns));

nv : number of values
xmin : start value of the list
xmax: end  value of the list

q:linspace(1,20,15);

q1:linspace(1,20,15),float;

xmin;
xmax;

q[length(q)],float;

delta (q)
=======
-----------------------------------------------------------
=>    Δq = float(q[i+1]-q[i]),  i=1,2,...,nv-1,nv
-----------------------------------------------------------

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

multiple plot : discrete data sets
========================
function :multiploti

website :
maxima:programming, different plots on the same canvas ,version 1.0
peter.vlasschaert@gmail.com,17/05/2020

1e) exact function pde : below
--------------------------------------------

f(x1,y1):=1/(2*π*x1)^(1/2)*exp(-y1^2/(2*x1));

w1:rhs(%);

rem : x1:time
          x2:space

2e) my multiple plot function.
---------------------------------------

a) sequence of values
===================
example : 0.1,0.2,0.3,........( step value:1)

multiploti(nu,g,xx,mx,tt,st):=block([m1:[]],
    xq:makelist(xx+mx,mx,0,10,0.1),
    tt:makelist(st*jj,jj,1,nu),  /*parameter*/
    for ii:1 thru nu do
                               (
        tq:tt[ii],
        yq:g(tq,xq),                              
        m1:endcons([discrete,xq,yq],m1)
                                 ),
        wxplot2d(m1)
)$

b) after some step value
===================
example : 0.1,0.3,.0.5....... ( step value :2)

makelist(ii,ii,1,10,2);  'list with step value 2'

multiplotis(g,xx,mx,ttt,stp):=block([m1:[]],
    xqq:makelist(xx+mx,mx,0,10,0.1),
      for ii:1 thru length(ttt) step stp do
                               (
        tqq:ttt[ii],
        yqq:g(tqq,xqq),                              
        m1:endcons([discrete,xqq,yqq],m1)
                                 ),
        wxplot2d(m1)
)$

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

numerical solution :pde, diffusion equation
================================

1e) pde : solve  
----------------------

r1:'diff(p,t,1) = 1/2*'diff(p,x,2);

BC :p(i,1) = 0, p(i,Nx) = 0                  →, i =all i ' time'
IC  : p(1,j) = gauss function = w1      →  j = all j 'space'
========================================
BC : (B)oundary (C)ondition
I  C : (I)nitial         (C)ondition
------------------------------------------------------------------------
------------------------------------------------------------------------
w1 : 'see above'
BC : see 7e)
------------------------------------------------------------------------

2e) pde :approximation ,numerical method : 'finite difference'
-----------------------------------------------------------------------------------
p = p[i,j]

r2:p[i+1,j]=1/2*'Δt/'(Δx)^2*(p[i,j+1]-2*p[i,j]+p[i,j-1])+p[i,j];

r2a:subst([i=i-1],r2);

r2a1:subst(['Δt =deltaT, 'Δx=deltaX],r2a);

scheme: 'finite difference'
------------------------------------

rem :'diff(p,t,1)  → forward difference in time
         'diff(p,x,2) → central difference   in space

3e) step :initial constant
----------------------------------

qq1:float(q[2]-q[1]);
qq2:float(q[3]-q[2]);
qq10:float(q[10]-q[9]);

Nt:1000;
Nx:100;
tmin:0.1;
tmax:10;
xmin : -5.0;
xmax:5.0;
t:linspace(tmin,tmax,Nt)$
x:linspace(xmin,xmax,Nx)$
deltaT:t[2]-t[1];
deltaX:x[2]-x[1];

rem :Δt = deltaT,Δx = deltaX

4e) step :create : zero matrix , Nt by Nx, with : Name 'p'
---------------------------------------------------------------------------
Menu → Algebra → Generate Matrix from Expression

info :p: genmatrix(lambda([i,j], 0), Nt, Nx);

p: genmatrix(lambda([i,j], 0), Nt, Nx)$

use matlab
========
p =zeros(Nt,Nx)

5e) analytical function solution : pde
-------------------------------------------------
nu : number of plots

t[1];

r2c:multiploti(1,f,x[1],Nx,t,t[1]);

use : matlab
==========
p(1,:) = analytical function(t(1),x) , *.m file ,' to start iteration scheme'.

t[1];

exact solution : pde
===============

r2d:multiplotis(f,x[1],Nx,t,45);

6e) p(1,j) = analytical solution ,start value
--------------------------------------------------------
→  firsrt values : for init time

for kk:1 thru  Nx do(
p[1,kk]:float(f(t[1],x[kk]))
)$

7e) iteration scheme
-----------------------------

for i:2 thru Nt do(    
    for j:1 thru Nx do(        
          if (  j=1 or j=Nx) then
          /*BC : Boundary Condition */
              p[i,j]:0
          else
         /*general iteration scheme */
             p[i,j]:(deltaT*(p[i-1,j+1]-2*p[i-1,j]+p[i-1,j-1]))/(2*deltaX^2)+p[i-1,j]         
                              )
                          ) $   
    

en[ii]:=makelist(p[ii,jj],jj,1,Nx);

first line matrix : p
=============

en[1]$

8e) plot 'matrix p'
------------------------

Matrix : multiplot: pp , step :stp
=======================

multiplotim(pp,xxm,ttm,stp):=block([m2:[]],
       em[ix]:=makelist(p[ix,jj],jj,1,Nx),
       for ii:1 thru length(ttm) step stp do
                               (                                   
        m2:endcons([discrete,xxm,em[ii]],m2)
                                 ),
        wxplot2d(m2)
)$

approximation solution : pde
======================

multiplotim(p,x,t,45);

numerical solution : 
==============

→ very difficult to check , phenomena : oscillation , blow up the best way : 
==================================================================
you need something to learning about ,"stability of a pde".
-----------------------------------------------------------------------------------------------------------------------
ref :Smith , Numerical Solution of Partial Differential Equations ,Oxford Press
-----------------------------------------------------------------------------------------------------------------------
→  pde :(p)artial (d)ifferential (e)quations

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