kill(all)$

biotech :appendix 2:asymptotic solution(second order diff. equation).
peter.vlasschaert@gmail.com,08/05/2020

ref : Fisher equation ( =  Nonlinear Diffusion-Reaction Phenomena)

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

load(lrats)$

option
=====
reason : rewrite 'numerics,floats' to 'rational numbers'
example : 0.5  = > 1/2

ratprint :false;

model :
=====

pp:ε*'diff(g,z,2)+'diff(g,z,1)+g*(1-g) =0;

BC :conditions :pg282
=================
---------------------------------------------------------------------------------
ref :MATHEMATICAL  BIOLOGY ,second,corrected edition
       → J.D. Murray
       → Springer VerLag
---------------------------------------------------------------------------------
---------------------------------------------------------------------------------
i=0,             g[0](-∞) = 1,g[0](∞) = 0,g[0](0) = 1/2           (BC1)
j=1,2,....      g[0](-∞) = 0,g[0](∞) = 0,g[j](0)  = 0      (i≠j)  (BC2)
---------------------------------------------------------------------------------

regular perturbation
===============

n:3$
su:0$
for i:0 thru n do  su: su+g[i](z)*ε^i$

order : peturbation series, n 
======================

su;

pp1:g(z)=su;

use : pp1 into 

p1:ε*diff(su,z,2)+diff(su,z,1)+su*(1-su);
p1:expand(p1);
for i:0 thru n  do(
eq[i]:coeff(p1,ε,i)=0
);

term : without ε
============
from p1

eq[0];

q1a:ode2(eq[0],g[0](z),z);

q1b:subst([z=0],q1a);

q1c:subst([g[0](0) =1/2],q1a);

q1d:logcontract(lhs(q1c))=rhs(q1c);

q1d1:rhs(q1d) = lhs(q1d);

BC1 : when  z=0, g[0](z) = g[0](0) = 1/2
==============================
find : %c

method  1: complex values
========
z = -1 = a+i*b = r*exp(i*θ)
=============================
r=sqrt((a)^2+(b)^2)=sqrt((-1)^2+(0)^2)
exp(i*π) =cos(π)+i*sin(π) = -1 +i*0
=============================
math :ln , maxima: log
=================
---------------------------------------------------
math       : ln(exp(i*π))       = i*π
maxima :  log(exp(%i*π)) = %i*π
---------------------------------------------------

r1a:subst([z=0,g[0](0)=1/2],q1d1);
r1a1:r1a,numer;

r1a2:subst([r1a1],q1d1);

ra1a3:exp(lhs(r1a2))=exp(rhs(r1a2));

ra1a4:ratsimp(solve(ra1a3,g[0](z))[1]);

method 2: I need real values
=======
need same form,find :a1 

rb1b:g[0](z) =(1-a1*exp(z))^(-1);

use : BC1

rb1b1:subst([z=0,g[0](0) = 1/2],rb1b);

find :a1

rb1b2:solve(rb1b1,a1);

rb1b3:subst(rb1b2,rb1b);

rb1b4:[rb1b3];

term : first order ε
==============
from : p1

eq[1];

q2a:args(lhs(eq[1]));

q2b:eq[1]-q2a[2];

use : info term : without ε  ,answer :method 2,see 'above'
------------------------------------------------------------------------------
use:q2b

q22a:subst(rb1b4,q2b);

q22b:ev(q22a,nouns);

q22c:lhs(q22b)-rhs(q22b) =0;

q22d:ode2(q22c,g[1](z),z);

=>  use model "pp",BC2

→ find , %c,from 'q22d'

q22d1:[z=0,g[1](0) =0];

q22d2:subst(q22d1,q22d);

q22d3:solve(q22d2,%c);

→subst, %c into 'q22d'

q22d4:subst(q22d3,q22d);

q22d4a:lhs(q22d4) =logcontract(rhs(q22d4));

q22d4b:lhs(q22d4) =factor(rhs(q22d4a));

term :second order ε^2
==================
from: p1

eq[2];

q3a:args(lhs(eq[2]));

q3b:eq[2]-q3a[2];

term :thrid order ε^3
==================
from:p1

eq[3];

q4a:args(lhs(eq[3]));
q4b:eq[3]-q4a[2];

asymptotic ( = solution)
==================
m*ε^2=O(ε^2),use 'pp1'
-----------------------------------------
→ g[2](z) = m
→ g[3](z) = 0

pp2a:makelist(g[i](z)=0,i,3,n);

pp2b:[g[2](z)=m];

pp2c:append(pp2a,pp2b);

pp2d:subst(pp2c,pp1);

use : rb1b3 & q224b

pp3e:[rb1b3,q22d4b];

pp4e:subst(pp3e,pp2d);