(%i1) | kill(all)$ |
ode : part7,intro: ricatti
[email protected],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
(%i2) |
load('contrib_ode)$ assume(t>0,CA>0,CB>0)$ |
equation (2) n=2,m=1, by 't=0 → CA=CA0'
(%i5) |
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); |
(%i7) |
eq4:subst(rhs(eq1),CA,eq3); eq5: lhs(eq3)=rhs(eq4); |
(%i8) | depends([y,f],[x]); $ |
(%i9) | eq6:'diff(y,x)+alpha*y=f; |
(%i10) |
depends([I],[x]); |
I is a integration factor.
(%i11) | eq7:ev('diff(I*y,x)=I*f,nouns,eq7); |
(%i12) | eq8:ev(eq7,I=%e^(alpha*x)); |
(%i13) | eq9:ev(eq8,nouns); |
make : eq5 = eq9
first part : eq5
x = t
y = CB
alpha = k2
(%i16) |
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
(%i18) |
eq13:factor(subst(rhs(eq5),f,eq12)); eq14:f=rhs(eq5); |
find integration factor I:eq5
(%i19) | dependencies; |
(%i20) | depends([I,f,CB],[t]); |
(%i21) | dependencies; |
(%i22) | 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
(%i23) | 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
(%i24) | eq17:ev(eq16,t=0,CB=0); |
(%i25) | eq18:solve(eq17,%c)[1]; |
(%i26) | eq19:ev(eq16,eq18); |
integration :integration(%e^(k2*t)*f,t),f see eq14
info : gamma_incomplete,see www
(%i27) | eq20:integrate(%e^(alpha*x)/x^2,x); |
make taylor expantion :
go to ' menu: Calculus → Get series.'
(%i28) | eq21:taylor(%e^(alpha*x)/x^2, x, 0, 4); |
(%i29) | eq22:integrate(eq21,x); |
Integration(I*f,t),use taylor expantion to do integration.
subst : tau=CA0*k1*t+1,alpha=k2
(%i30) | eq23:tau=part(eq14,2,2,1); |
(%i31) | eq24:diff(rhs(eq23),t); |
(%i32) | eq25:k2*solve(eq23,t)[1]; |
build %e^(k2*t) ⇒ function of tau
(%i33) | eq26:%e^(lhs(eq25))=%e^(rhs(eq25)); |
build : %e^(k2*t)*f, f=f(tau)
(%i34) | eq27:ratsubst(tau,rhs(eq23),eq14); |
(%i35) | eq28:v=integrate(rhs(eq25)*rhs(eq27),tau); |
(%i37) |
eq29:part(eq28,2,3,1); eq30:ratsubst(rhs(eq23),lhs(eq23),eq29); |
(%i38) | taylor(eq30, t, 0,3); |
(%i39) | eq31:CB=1/rhs(eq15)*rhs(eq28)+%c/rhs(eq15); |
t = 0 ,tau = 1+CA0*k1*t,tau = 1 and CB=0
(%i45) |
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)); |