kill(all)$

module 3a: solve system : differential equations. 2 by 2
solve : coolinglaw newton (laplace transform)
peter.vlasschaert@gmail.com , 11/10/2017

assume (k>0)$

def:cooling law :Newton.    
      
   'diff(T,t,1)=-(T-Te)*k

T0 = initial temperature ( initial state ),see below
Te = ambient temperature ( final   state ),see below
T   = temperature ' function of time ' = T( time ) = T(t);
rem : IC => T0 (Initial Condition)

L{T(t)} = Ts  =>    T(t) = L(-1)L{T(t)} = L(-1){Ts}
L        = Laplace Transform
L(-1)  = Inverse Laplace Transform

depends(T,t);

pT1:diff(T(t),t)=-k*(T(t)-Te);
pTic : T(0) =T0;

first :  Laplace Transform  'pT1'

a)  L{'diff(T,t,1)} =  s*Ts -T(0)

p1:Ts=laplace(T(t), t, s);
p2:laplace(diff(T(t),t), t, s);
p3:ratsubst(lhs(p1), rhs(p1), p2);

b)  L{T(t)} = Ts ,use  'a)'  for 'pT1'  : '  L transform '
rem : L{const} =const/s

pT2:laplace(pT1,t,s);
pT2a : ev(pT2,rhs(p1)=lhs(p1));

c) solve for  Ts = ?
menu:Equations:Solve.

p4:solve(pT2a,Ts)[1];

d) solution : L(-1){Ts} = T(t)

p5:T(t)=ilt((T(0)*s+Te*k)/(s^2+k*s), s, t);

find : ic ,t=0

q1:'limit(part(p5,1),t,0)='limit(part(p5,2),t,0);
q2:'limit(part(p5,1),t,0)=limit(part(p5,2),t,0);
q3:ev(q2,pTic);

=> T0 initial temperature

find : steady state , t=infinity  (time) 

q11:'limit(part(p5,1),t,inf)='limit(part(p5,2),t,inf);
q22:'limit(part(p5,1),t,inf)=limit(part(p5,2),t,inf);

after a long time  takes same temperature = Te = ambient temperature

check solution : 'p5' solve 'pT1'

p6:ev(pT1,p5,pTic,diff);
p7:is(p6);

kill(all);
depends(T,t);

three zones cooling : Te = put in egg in boiling water = 'Te' of 'Te=100*e^(-kt)'
                                         T2 = egg white ( = Albumen)      = T[1]   at (t=0) => 40
                                         T1 = yolk                                         = T[2]   at (t=0) => 45

w1:Te =0;   /*ice bad */
w2:a=2;      /*no physical meaning  a,b */
w3:b=3;

ps[1]:'diff(T[1],t)=a*(T[2]-T[1]);
ps[2]:'diff(T[2],t)=a*(T[1]-T[2])+b*(Te-T[2]);

pss[1]:expand(ev(ps[1]));
pss[2]:expand(ev(ps[2]));

build coefficient  2 by 2 :ODE

M1:zeromatrix(2,2)$
for i:1 thru 2 do
    for j:1 thru 2 do
  M1[i,j]:coeff(part(pss[i],2),T[j]);;

coefficient matrix ps[1],ps[2]

M1;

T[1] =>  f(t)
T[2] =>  g(t)

depends([f],t)$
depends([g],t)$

w4:diff(f(t),t) = M1[1,1]*f(t)+M1[1,2]*g(t);
w5:diff(g(t),t) = M1[2,1]*f(t)+M1[2,2]*g(t);

find : L from 'w4 and w5'

w4s:laplace(w4,t,s)$
w5s:laplace(w5,t,s)$
w6:laplace(f(t),t,s) =fs$
w7:laplace(g(t),t,s)=gs$
w6a:ev(w4s,w6,w7);
w7a:ev(w5s,w6,w7);

solve :  {fs,gs}

w8a:solve([w6a,w7a],[fs,gs])[1][1];
w8b:solve([w6a,w7a],[fs,gs])[1][2];

find : L(-1) get   {f(t),g(t)}

w9a:ev(part(w8a,2),w1,w2,w3)$
w9b:ev(part(w8b,2),w1,w2,w3)$
w10a:f(t)=ilt(w9a,s,t);
w10b:g(t)=ilt(w9b,s,t);

insert : initial conditions 'see above' (t=0,time),{f(t),g(t)}
f(0)  = 40
g(0) = 45

wx1:f(0) = 40;
wx2:g(0) = 45;
wa:ev(w10a,wx1,wx2);
wb:ev(w10b,wx1,wx2);

check : ode 2 by 2, {w4,w5}

waa:ratsimp(ev(w4,wa,wb,diff,w1,w2));
wbb:ratsimp(ev(w4,wa,wb,diff,w1,w2));
is(waa);
is(wbb);