kill(all)$
load(eigen)$

module 3:cooling law 'system differential equations'. 
peter.vlasschaert@gmail.com ,           25/03/2017

cooling law :Newton.
T0 = initial temperature ( initial state ),see below
Te = ambient temperature ( final   state ),see below
T  = temperature ' function of time ' = T( time ) = T(t)

assume(k>0);

pT1:'diff(T,t)=-k*(T-Te);

p1:ode2(pT1,T,t);

How to get '%c' => t=0 -> T=T0

p2:ev(p1,t=0);
p3:solve(p2,%c)[1];
p4:ev(p3,T=T0);

p5:ev(p1,p4);

check : of the initial state of the cooling process.

find initial state 'p5' ' t-> 0 ',time domain
( Laplace transformation , s domain )

p6:limit(part(p5,2), t, 0);

p7:T=ratsimp(part(p5,2));

find steady state 'p7' ' t-> infinity ',time domain
( Laplace transformation , s domain )

p8:limit(p7,t,inf);

kill(all);

three zones cooling : Te = put in egg in boiling water = 'Te' of 'Te=100*e^(-kt)'
                      T2 = egg white ( = Albumen)      = T[1]
                      T1 = yolk                        = T[2]
-----------------------------------------------------------------------------------------------------
http://www.ltz.de/en/news/lohmann-information/An-egg-a-day-the-physiology-of-egg-formation.php ,'fig'
-----------------------------------------------------------------------------------------------------
temperature range : Te(t) or constant , T2(t) , T1(t)
-----------------------------------------------------------------------------------------------
egg model :
-----------------------------------------------------------------------------------------------
mitopencourseware : 
Instructor/speaker: Prof. Arthur Mattuck
Home » Courses » Mathematics » Differential Equations » Video Lectures » 
-----------------------------------------------------------------------------------------------
-----------------------------------------------------------------------------------------------
https://ocw.mit.edu/courses/mathematics/18-03-differential-equations-spring-2010/video-lectures
/lecture-24-introduction-to-first-order-systems-of-odes/
-----------------------------------------------------------------------------------------------

a,b : no physical meaning, 'different conductions = different regios'
Te = 0    (ice bad,course mit use this)
Te = 100  (boiling egg)

edit : Te:0  (course mit) -> Te:100

Te :0;

depends(T,t);

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

a:2;
b:3;

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

eliminate : T2 

p9:solve(pT1,T[2])[1];
a:2;
b:3;

p10:ratsimp(ev(pT2,p9));

p11:ratsimp(ev(lhs(p10)-rhs(p10)=0,nouns));

find : T
        1

rem : characsteristic eq: r^2+7*r+6=0,roots = -1,6

p12:ode2(p11,T[1],t);

intial conditions T[1](t=0)=40
                  T[2](t=0)=45 

find : T
        2

p13:ev(p12,%k1=c1,%k2=c2);

T[1](t=0)=40,initial temperature (yolk)

p131:ev(p13,T[1]=40,t=0);

p14:ev(p9,p13);

p15:ev(p14,nouns);

ratsimp(ev(p14,nouns));

T[2](t=0)=45,initial temperature (egg white)

p151:ev(p15,T[2]=45,t=0);

how to do this :menu
Equations -> solve linear system -> ask questions
 2 -> first   field : p131
      second  field : p151
      third   field : c1,c2

p16:linsolve([p131, p151], [c1,c2]);

p17:p16[1];
p18:p16[2];

T (t) = temperature 'yolk with time'.
 1

p19:ev(p13,p17,p18);

T (t) = temperature 'white of the egg with time'.
 2

p20:ratsimp(ev(p15,p17,p18));

kill(all);

solve : the above problem again ,with eigenvectors
                                 with eigenvalues.
general method :solving system of differential equations.
X = T[1]
Y = T[2]
mit course:
https://ocw.mit.edu/courses/mathematics/18-03-differential-equations-spring-2010/video-lectures
/lecture-25-homogeneous-linear-systems-with-constant-coefficients/

depends(X,t);
depends(Y,t);

pd1:'diff(X,t) = -2*X+2*Y;
pd2:'diff(Y,t) = 2*X-5*Y;

what with : general 'solution'
lambda = parameter
X= a1*e^(lambda*t)
Y= a2*e^(lambda*t)

p21:X=a1*%e^(lambda*t);
p22:Y=a2*%e^(lambda*t);

p23:ev(pd1,p21,p22);
p24:ev(pd2,p21,p22);

p231:ev(p23,nouns)/%e^(t*lambda);
p241:ev(p24,nouns)/%e^(t*lambda);

p25:ratsimp(p231);
p26:ratsimp(p241);

p27:rhs(p25)-lhs(p25)=0;
p28:rhs(p26)-lhs(p26)=0;

p29:coeff(part(p27,1),a1);
p30:coeff(part(p27,1),a2);
p31:coeff(part(p28,1),a1);
p32:coeff(part(p28,1),a2);

build matrix M1 ' coefficient matrix '

M1:zeromatrix(2,2);

M1[1,1]:p29$
M1[1,2]:p30$
M1[2,1]:p31$
M1[2,2]:p32$

coefficient matrix,p27,p28

M1;

call :characteristics polynomial of a matrix
module 1 : charpoly form ' load(eigen)'

p33:ratsimp(determinant(M1))=0;

lambda-sum*lambda+product  = 0 'quadratic equation '
sum     = lambda1+lambda2
product = lambda1*lambda2
------------------------------
example: lambda1=-6,lambda2=-1

p34:solve(p33,lambda);

p35:p34[1];
p36:p34[2];

p351:part(p35,2);
p361:part(p36,2);

find eigenvector corresponding to 'lambda=-1'

p37:a1*M1[1,1]+a2*M1[1,2]=0;
p38:a1*M1[1,1]+a2*M1[1,2]=0;

subst : lambda = -1 'p36' in to p37,p38

p39:ev(p37,p36);
p40:ev(p38,p36);

menu :equations ->solve linear system 'p39,p40' for a1,a2

eigenvectors (  lambda = -1)
module 1 : eigenvectors ,uniteigenvectors

p41:linsolve([p39, p40], [a1,a2]);

p42:ev(p41,%r1=1);

p43:p42[1];
p44:p42[2];

subst : lambda = -6 'p35' in to p37,p38

p45:ev(p37,p35);
p46:ev(p38,p35);

menu :equations ->solve linear system 'p45,p46' for a1,a2

eigenvectors (  lambda = -6)
module 1 : eigenvectors ,uniteigenvectors

p47:linsolve([p45, p46], [a1,a2]);

p48:ev(p47,%r2=-2);

p49:p48[1];
p50:p48[2];

superposition principle :C1,C2
before : C11 = c1/2,C22 = c2

p51:X=C11*part(p43,2)*%e^(p361*t)+C22*part(p49,2)*%e^(p351*t);
p52:Y=C11*part(p44,2)*%e^(p361*t)+C22*part(p50,2)*%e^(p351*t);

general system : without values. 'appendix'