kill(all)$

ode : part7,intro: ricatti 
peter.vlasschaert@gmail.com,06/02/2017

# Theory Ricatti : ODE (reaction kinetics)

(dy/dx) = p(x)*y^2+q(x)*y+r(x)
Why: Consecutive Reactions and Ricatti
                      k1   k2
reaction kinetics : A -> B -> C
dCA/dt = - k1*CA^n          (1)
dCB/dt = k1*CA^n-k2*CB^m    (2)
dCC/dt = k2*CC^l            (3)
                                              
integer values :  n,m,l 
concentration  : CA,CB,CC

load('contrib_ode)$
assume(t>0,CA>0,CB>0)$

equation (2) n=2,m=1, by 't=0 -> CA=CA0'

eq1: CA = CA0/(1+k1*CA0*t);
eq2:'diff(CB,t)=k1*CA^2-k2*CB^1;
eq3:lhs(eq2)-part(eq2,2,2)=rhs(eq2)-part(eq2,2,2);

eq4:subst(rhs(eq1),CA,eq3);
eq5: lhs(eq3)=rhs(eq4);

depends([y,f],[x]);                                                                        $

eq6:'diff(y,x)+alpha*y=f;

depends([I],[x]);

I is a integration factor.

eq7:ev('diff(I*y,x)=I*f,nouns,eq7);

eq8:ev(eq7,I=%e^(alpha*x));

eq9:ev(eq8,nouns);

make : eq5 = eq9

first part : eq5

x = t
y = CB
alpha = k2

eq10:subst(t,x,eq9);
eq11:subst(CB,y,eq10);
eq12:subst(k2,alpha,eq11);

second part :eq5

f=(CA0^2*k1)/(CA0*k1*t+1)^2

eq13:factor(subst(rhs(eq5),f,eq12));
eq14:f=rhs(eq5);

find integration factor I:eq5

dependencies;

depends([I,f,CB],[t]);

dependencies;

eq15:I=part(eq13,1,2);

seperated,find CB=CB(t) from eq5
solve eq5:
I = I(t),f=f(t),%c = integration constant
formule : CB = 1/I*integration(I*f,t)+%c/I

eq16:CB =ev(1/rhs(eq15)*integrate(rhs(eq15)*rhs(eq14),t)+%c/rhs(eq15),nouns);

find %c , same part6, t=0 and CB=0 in eq16

eq17:ev(eq16,t=0,CB=0);

eq18:solve(eq17,%c)[1];

eq19:ev(eq16,eq18);

integration :integration(%e^(k2*t)*f,t),f  see eq14

info : gamma_incomplete,see www

eq20:integrate(%e^(alpha*x)/x^2,x);

make taylor expantion :
go to  ' menu: Calculus -> Get series.'

eq21:taylor(%e^(alpha*x)/x^2, x, 0, 4);

eq22:integrate(eq21,x);

Integration(I*f,t),use taylor expantion to do integration.

subst :  tau=CA0*k1*t+1,alpha=k2

eq23:tau=part(eq14,2,2,1);

eq24:diff(rhs(eq23),t);

eq25:k2*solve(eq23,t)[1];

build %e^(k2*t) => function of tau

eq26:%e^(lhs(eq25))=%e^(rhs(eq25));

build : %e^(k2*t)*f, f=f(tau)

eq27:ratsubst(tau,rhs(eq23),eq14);

eq28:v=integrate(rhs(eq25)*rhs(eq27),tau);

eq29:part(eq28,2,3,1);
eq30:ratsubst(rhs(eq23),lhs(eq23),eq29);

taylor(eq30, t, 0,3);

eq31:CB=1/rhs(eq15)*rhs(eq28)+%c/rhs(eq15);

t = 0 ,tau = 1+CA0*k1*t,tau = 1 and CB=0

eq32:ratsubst(lhs(eq23),rhs(eq23),eq31);
eq33:solve(eq23,t)[1];
eq34:ev(eq32,tau=1,CB=0,lhs(eq33)=rhs(eq33));
eq35:solve(eq34,%c)[1];
eq36:ev(eq32,eq33,eq35);
eq37:factor(ev(eq36,eq23));