kill(all)$

load packages :

load(lrats)$

heat : nonhomogeneous problem (bc (t)≠ 0),(ic = f(x),@t=0),version 1.0 (part 3)
peter.vlasschaert@gmail.com,25/07/2019.

rem : no heat generation  ( = without a source)

ref : elementary applied partial differential equations
        with fourier series and boundary value problems
        (richard haberman,prentice hall)

Model  3 ( cartesian coordinates : u(x,t), unsteady = du/dt ≠ 0)
=======
x: space 
t : time         > 0
u : physical quantity ( here : temperature)
D : constant  > 0
l   : length ( for example : rod) , 0<x<l

conditions
---------------

pde :D* 'diff(u,x,2) = 'diff(u,t) 
bc :  u(0,t)  =A(t) time dependent
        u(l,t  ) =B(t)  time dependent
ic  :  u(x,0) = f(x)

type : linear homogeneous partial differential equation.
           bc : linear nonhomogeneous boundary conditions.

How to solve  => separation of variables.

Equilibrium temperature : uE
===================

assume(D>0);

assume(l>0);

p1:D*diff(u(x,t),x,2) = diff(u(x,t),t);

independent of time only function of space
----------------------------------------------------------

depends(uE,x);

dependencies;

u(x,t) => uE(x)
===========

1e) independent of time: du/dt =0
----------------------------------------------

p2:lhs(p1) = 0;
p2a:lratsubst([D=1],p2);
p2b:lratsubst([u(x,t)=uE(x)],p2a);

2e) How to solve uE for bc
-------------------------------------

p3:desolve(p2b,uE(x));

coef : x
====

p3a: coeff(rhs(p3),x)=Cx;

p3b:lratsubst(p3a,p3);

coef :constant
====

p3c:ev(rhs(p3b),x=0)=C;

p3d:lratsubst([p3c],p3b);

3e) how to find time dependent
-------------------------------------------
uE(0) = A(t)
uE(l)   =B(t)

p4a:lratsubst([x=0],p3d);
p4b:lratsubst([x=l],p3d);

4e) solve for Cx and C
--------------------------------

p5:solve([p4a,p4b],[Cx,C]);

p5a:flatten(p5);

5e) equilibrium temperature distribution
-------------------------------------------------------

p6:lratsubst(p5a,p3d);

p6a:[uE(0) = A(t),
uE(l)   = B(t)];

p6b:lratsubst(p6a,p6);

p6c:expand(p6b);

rem : uE(x)  => uE(x,t)
================

p6d:lratsubst([uE(x)=uE(x,t)],p6c);

dependencies;

general relation : u(x,t) = v(x,t)+uE(x)
============================
v(x,t) : solution to problem with homogeneous bc conditions

model v(x,t)
=========
pde :D* 'diff(v,x,2) = 'diff(v,t) + ?
bc :  v(0,t) =0    (homogeneous bc)
        v(l,t  ) =0    (homogeneous bc)
ic  :  v(x,0)   ?

see : part 1  , u(x,t) =>v(x,t)
====================

p7:lratsubst([u(x,t)=v(x,t)],p1);

p7a:u(x,t) = v(x,t)+uE(x,t);

p7b:expand(lratsubst([p6d],p7a));

subst : 'p7b' into p1

p7c:lratsubst([p7b],p1);

p7d:ev(p7c,nouns);

p7e:ψ(x,t) =rhs(p7d-diff(v(x,t),t));

meaning : v(x,t) : displacement from quilibrium
===================================

p7f:solve(p7a,v(x,t))[1];

solution v(x,t)
==========
see :part 1:use p7c

p8:v(x,t)=sum(C4[n]*%e^(-(%pi^2*D*n^2*t)/l^2)*sin((%pi*n*x)/l),n,1,∞);

Why :see: part 2,website(also)
------------------------------------------

p8a:args(rhs(p8))[1];

find coef :

p8b:ratcoeff(p8a,sin((%pi*n*x)/l)) = C5(t);

we can rewrite p8a, use orthogonality relation for sine function  => C5(t)
--------------------------------------------------------------------------------------------------

p8c:lhs(p8) =intosum(sum(C5(t)*sin((%pi*n*x)/l),n,1,∞));

ic : condition (@t=0)
================

p9:lratsubst([t=0],p7a);

ic: from model 2: u(x,0) = f(x)
======================

p10a:lratsubst([u(x,0)=f(x)],p9);

p10b:solve(p10a,v(x,0))[1];

find : uE(x,0)
----------------

p10c:lratsubst([t=0],p6d);

find : v(x,0),ψ(x,t)
-----------------------

p10d:subst([p10c],p10b);

p10e:p7e;

model v(x,t)
=========
pde :D* 'diff(v,x,2) = 'diff(v,t) + ψ(x,t)
bc :  v(0,t) =0    (homogeneous bc)
        v(l,t  ) =0    (homogeneous bc)
ic  :  v(x,0) , see above
        ψ(x,t)  ,see above

p10f:D*diff(v(x,t),x,2) = diff(v(x,t),t) + ψ(x,t);

p10g:intosum(lratsubst([p8c],p10f));

evaluation of: 'above
------------------

p10h:intosum(ev(p10g,nouns));

simplify :p10h

solution : p10h
===========

1e) ψ(x,t) :odd extend ,period  2*L, 'periodic function',expand fourier series

fourier series
===========

p11a:lhs(p7e) = intosum(sum(ψn(x,t)*sin((n*π*x)/l),n,1,∞));

find coef:sin((%pi*n*x)/l)

p11b:args(rhs(p11a));

p11c1:part(p11b[1],1);

p11c2:part(p11b[1],2);

fourier coef.
=========

p11c:p11c1=2/l*'integrate(lhs(p11a)*p11c2,x,0,l);

2e) input  ψ(x,t) into 'p10h'
-----------------------------------

p11d:lratsubst([p11a],p10h);

p11d1:lhs(p11d)-rhs(p11d) =0;

p11d2:intosum(p11d1);

p11d3:sumcontract(p11d2);

p11d4:factor(p11d3);

3e)  How to find coef : sin((%pi*n*x)/l) ≠ 0
--------------------------------------------------------

p11e:lhs(p11d4);

p11e1:part(p11e,1,1,1,1)=0;

p11e2:solve(p11e1,ψn(x,t))[1];

4e)  How to find coef : c5(t)
--------------------------------------------------------

p11e3:ode2(p11e2,C5(t),t);

copy and paste: output 'p11e3',integration 0 to mm =t , maxima 'problem'

p11e3a:C5(t)=%e^(-(%pi^2*D*n^2*t)/l^2)*(%c-'integrate(%e^((%pi^2*D*n^2*t)/l^2)*ψn(x,t),t,0,mm));

final result :find u(x,t)
================

find : v(x,t)
========
use : p8c

p12a:v(x,t) = lratsubst([p11e3a],rhs(p8c));

p12b:expand(p12a);

How to find  : %c  => from ic : condition @t=0
--------------------------------------------------------------

a) reason mm=0
------------------------

p12c:lratsubst([mm=0],rhs(p12b));

b) use t=0 into p12c
---------------------------

p12c1:intosum(lratsubst([t=0],p12c));

c) orthogonality,v(x,0)
------------------------------
use :p10d

p12c2:rhs(p10d) = p12c1;

d) fourier coef: %c
-----------------------------

p12c3:%c = 2/l*integrate(rhs(p10d)*p11c2,x,0,l);

=================
final final result : u(x,t)
=================

  =====================================================================

p13:p7b;
p13a:lratsubst([mm=t],p12b);
p13b:p12c3;

=====================================================================

rem : integral ( same value : n=1,2,3,....) for orthogonality
==============================

p14c1:'integrate((sin((%pi*n*x)/l))^2,x,0,l) =lratsubst([n=1],integrate((sin((%pi*n*x)/l))^2,x,0,l));

p14c2:'integrate((sin((%pi*n*x)/l))^2,x,0,l) =lratsubst([n=5],integrate((sin((%pi*n*x)/l))^2,x,0,l));

1e) integral useful:energy or power 
2e) orthogonality of sines