kill(all)$

biotech: method of lines 'fisher nonlinear pde', version 1
peter.vlasschaert@gmail.com,26/05/2020

see:website:
biotech:apppendix : method of lines 'fisher nonlinear pde', version 1
peter.vlasschaert@gmail.com,26/05/2020

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

load(lrats)$

ratprint:false;

→  method of lines (MOL method):need 'system of discrete odes'
==================================================
'see below'

model :fisher nonlinear pde
======================

q1:diff(u(x,t),t) = D*diff(u(x,t),x,2)+r*u(x,t)*(1-u(x,t));

D= diffusion coef.
r  = reaction rate
u = u(x,t)
---------------------------------------------------------------------------
→ : solve 'q1' ,need BC,IC
----------------------------------------------------------------------------
BC : Boundary Condition
IC   : Initial         Condition
-----------------------------------------------------------------------------
domain

 ===========================================

info 1:EqWorld
========================
Exact Solutions > Nonlinear Partial Differential Equations
Fisher equation:'diff(u(x,t),t,1)=D*('diff(u(x,t),x,2))+r*(1-u(x,t))*u(x,t), D=1,C=1
sol1 :w(x,t)=1/[1 +C*exp(−5/6*r*t±1/6*√(6*r)*x)]^2
---------------------------------------------------------------------------------
http://eqworld.ipmnet.ru,first solution.(+ solution),we will use.

def :f(x,t,r,C) 
==========

f(x,t,r,C) :=1/(1 +C*exp(−5/6*r*t+1/6*sqrt(6*r)*x))^2;

→ BC    0 ≤ t≤ tfin
----------------------------

p1a:u(0,t) =f(0,t,1,1);
p1b:u(1,t) =f(1,t,1,1);

rem : tfin = time for simulation the problem.

→ iC(t=0)    0 ≤ x≤ L =1  
--------------------------------

p1c:u(x,0) = f(x,0,1,1);

(m)ethod (o)f (l)ines(=MOL) ,( nonlinear pde)
==================================

→ first step : discrete version for ,'q1'
---------------------------------------------------

eq[i]:='diff(u[i],t) = D/h^2*(u[i+1]-2*u[i]+u[i-1]+r*u[i]*(1-u[i]));

info 2 :rk4
=======
Function: rk  
    rk (ODE, var, initial, domain)  
    rk ([ODE1, …, ODEm], [v1, …, vm], [init1, …, initm], domain)
example :
=======
→  dx/dt = 4-x^2-4*y^2  (1)   
      dy/dt = y^2-x^2+1     (2)
For t between 0 and 4, and with values of -1.25 and 0.75 for x and y at t=0: 
→  sol: rk([4-x^2-4*y^2,y^2-x^2+1],[x,y],[-1.25,0.75],[t,0,4,0.02])$
The plot will show the solution for variable y as a function of t. 
→   wxplot2d ([discrete,makelist([p[1],p[3]],p,sol)], [xlabel,"t"],[ylabel,"y"])$

→ second step : find ODEm, m =10
--------------------------------------------------

p2:makelist(eq[ii],ii,1,10);

p2a:makelist(rhs(p2[jj]),jj,1,length(p2));

→ thrid step : 0 ≤ xx≤ L =1  &  parameters
--------------------------------------------------------
find :h ?

L:1; /*dimension domain x-direction*/

p3a:makelist(L*kk/(length(p2)+1),kk,0,length(p2)+1)$

xx:p3a;

→ parameters

h:xx[2]-xx[1];
r:1;
D:1;

p3b:ev(p2a,nouns);

→fourth step : list of variable:
---------------------------------------

p4a:listofvars(p3b);

BC : u[0],u[11],different for time step :Δt, found fom BC, 'above'

p4b:makelist(p4a[i],i,2,length(p4a)-1);

→fifth step :u[0],u[11],different for time step :Δt
----------------------------------------------------------------
tbegin=0
tend   =tfinal

u[0]:rhs(p1a);
u[length(p4a)-1]:rhs(p1b);

p4c:ev(p3b,nouns);

→sixth step :IC condition (when t=0)
----------------------------------------------------------------
use: x,p1c

p5a:xx;
p5b:makelist(lratsubst([x=xx[ll]],rhs(p1c)),ll,1,length(p5a));

p5bb:float(p5b);

1e) rk :not allowed is index variables
2e) rk : differential equations mayt not contain indepent variable (=time)
--------------------------------------------------------------------------------------------------
see : 'use rk:info 2', can't be used.

p6:makelist(u[i] = concat(u,i),i,1,length(p2));

p6a:lratsubst(p6,p4c);
p6b:lratsubst(p6,p4b);

v:50;   /*number of time steps */

pp:makelist(makelist(0,j,1,v),i,1,10)$

ui  → pp[i] list  ,i =1...number of  differential equations = 10

try compute with : eulermethod

fpprintprec : 3;

ratio=L*Δt/h;

matrix : a
=======

a:transpose(apply(matrix,pp))$

→  IC , will be first row fill with values p5bb

p7:makelist(a[1,i]:p5bb[i],i,1,length(p2))$
a$

p7a:makelist(p6b[i]=h*i,i,1,length(p2));

result :system of differential equations
=============================

p7b:p6a,float;

→  time :step

Δt:0.002;

→ Courant number =ratio
=====================
       L   =length domain
       h   = space division
       Δt = time step

ratio=L*Δt/h;

use : p7b , t →i*Δt
===============

p7bb:lratsubst([t=i*Δt],p7b);

implementation :euler method
=======================

for i:1 thru v-1 do(
    d1:float(makelist(concat(u,ii)=a[i,ii],ii,1,length(p2))),
   for j: 1 thru length(p2) do
    (               
      d2: float(lratsubst(d1,p7bb[j])),
         a[i+1,j]:a[i,j]+Δt*ev(d2 ,nouns)
        )                                              
);

result : euler algo 'above'
=====::=============

a:float(a);

plot : matrix 'a', u[i,j]
===============
i→x : horizontal                                  ( i = time evolution (i*Δt))
j→u : discrete values for fisher pde ( j = space division        )

wxdraw2d(image(a,0,0,length(p2),v),dimensions=[1500,1500]);

wxdraw(gr2d(colorbox =false, image(a,0,0,v,length(p2)))
                    ,dimensions = [1000,1000]  
                    );

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