kill(all)$

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

load packages
===========

load(lrats);

assume(D>0);

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

steady state
==========

p2:lhs(p1)=0;

with general : Q(x)

p2a:ode2(p2,u,x);

p2b:subst([Q(x) =m1*x+m2],p2);

example :Q(x) = m1*x+m2

1e)  p2b
=======

p2b1:ode2(p2b,u,x);

2e) p2a
======

p2a1:subst([Q(x) =m1*x+m2],p2a);

p2a2:ev(p2a1,nouns);

p2a3:expand(p2a2);

Model  4 ( 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) + Q(x)
bc :  u(0,t)  =A 
        u(l,t  ) = B  
ic  :  u(x,0) = f(x)

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

rem:: steady state solution  = v(x)     , t → ∞
          transient      soluttion  = w(x,t)  , t

1e) general solution
----------------------------

p3:u(x,t) = v(x)+w(x,t);

2e) 'p3' into 'p1'
------------;--------

a)  u → u(x,t)
------------------

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

b)  subst
------------

p4a:lratsubst([p3],p4);

p4b:ev(p4a,nouns);

p4c:expand(p4b);

3e) use bc : (@ x=0,x=l)
---------------------------------

p5a:lratsubst([x=0],p3);
p5b::lratsubst([x=l],p3);

p5a1:A=rhs(p5a);
p5b1:B=rhs(p5b);

4e)  use ic :(@t=0)
--------------------------

p6a:lratsubst([t=0],p3);

p6b:f(x) =rhs(p6a);

steady state : solution
================
w(x,t) = 0  => 'p4c'

p7:lratsubst([w(x,t)=0],p4c);

p7a:lratsubst([v(x)=v],p7);

p7b:ode2(p7a,v,x);

how to find : %k1 and %k2
-------------------------------------
use : set :v(0) = A, v(l) = B

p8:[v(0)=A,v(l) = B];

p8aa:p8[1];
p8bb:p8[2];

use :p5a1,p5b1

@ x=0
---------

p8a:subst(p8,p5a1);

p8a1:solve(p8a,w(0,t));

@ x=l
---------

p8b:subst(p8,p5b1);

p8b1:solve(p8b,w(l,t));

use : Q(x) =m1*x+m2
------------------------------
use : p7b,p8

p9:[Q(x)=m1*x+m2];

p9a:subst(p9,p7b);

p9b:expand(ev(p9a,nouns));

use : bc from , 'p8'
-------------------------

p10a:rhs(p8aa) =subst([x=0],rhs(p9b));
p10b:rhs(p8bb) =subst([x=l],rhs(p9b));

result : %k1,%k2
=============

p10c:solve([p10a,p10b],[%k1,%k2]);

p10c1:flatten(p10c);

result : v(x)
========;

p10c2:subst(p10c1,p9b);

model : w(x,t)
------------------
bc: p8a1,p8b1
ic: p6b

pde :D* 'diff(w,x,2) = 'diff(w,t)   : 0<x<l
bc :  w(0,t)  =0 
        w(l,t  ) = 0  
ic  :  w(x,0) =?

pde: p4c,p7a
-------------------

p11:lratsubst([v(x)=0,Q(x)=0],p4c);

bc,ic

p11a:p8a1;
p11b:p8b1;
p11c:p6b;

p11c1:solve(p11c,w(x,0));

see: part1,part2,part3

general solution: w(x,t)
==================
separation variables 

p12:w(x,t) =sum(c(n)*%e^(-(%pi^2*D*n^2*t)/l^2)*sin((%pi*n*x)/l),n,1,inf);

p12a:subst([t=0],p12);

use :p11c1

p12b:lratsubst(p11c1,p12a);

assume(l>0);

orthogonality : sine, fourier series sine
-----------------------------------------------------

p12c:c(n)=2/l*'integrate(lhs(p12b)*sin((%pi*n*x)/l),x,0,l);

final solution : u(x,t),Q(x) = m1*x+m2
============================

p13a:p3;

p13b:v(x) =rhs(p10c2);

p13c:p12;

p13d:p12c;

p13:subst([p13b,p13c,p13d],p13a);

rem:bc=constant  => 'see above' ,u(x,t)=w(x,t)+v(x)
        bc = f(t)          =>                      ,u(x,t)=w(x,t)+v(x,t)
=======================================
bc: boundary conditions
ic :  initial       conditions
=======================================