Maple Code

Example 11.1-2. Heating of a Finite Slab

> restart;

Equation for heat conduction in a solid, after the insertion of Fourier's Law. (Eq. 11.1-2 for heat conduction in y direction)

> eqn1:=diff(T(y,t),t)-alpha*diff(T(y,t),y$2)=0;

Switch to dimensionless variables with dchange.

> with(PDEtools,dchange):

The first transformation is from Eq. 11.1-11. It is written so that Theta will correspond to the initial condition and two boundary conditions given in 11.1-15 and 11.1-16.

> transform:={T=T1-Theta*(T1-T0),y=b*eta,t=b^2*tau/alpha};

> eqn2:=dchange(transform,eqn1,[Theta(tau,eta),tau,eta]);

Eq. 11.1-14

> simplify(-eqn2*b^2/alpha/(T1-T0));

As Maple shows us, this equation can be solved by separation of variables because the solution depends on two linearly independent functions of one variable each.

> pdsolve(eqn2);

Defining Theta(eta,tau) as a product of the two linearly independent functions,

> Theta:=(eta,tau)->f(eta)*g(tau);

> eqn3:=diff(Theta(eta,tau),tau)=diff(Theta(eta,tau),eta$2);

> eqn4:=simplify(eqn3/(f(eta)*g(tau)));

Eqn4 is only valid if both sides equal a constant, which we will call -c^2.

> left:=lhs(eqn4)=-c^2;

> right:=rhs(eqn4)=-c^2;

We now solve each independently.

> gtau:=rhs(dsolve(left,g(tau)));

> feta:=rhs(dsolve(right,f(eta)));

In order to be consistent with the book, we will rename the constant of integration for gtau and feta.

> gtau:=subs(_C1=A,gtau);

> g:=unapply(gtau,tau);

> feta:=subs({_C1=C,_C2=B},feta);

> f:=unapply(feta,eta);

By symmetry, everything that goes on at the ends of the slab must be equal. This is shown in B.C.'s 1 and 2, since both are zero.

> eqn5:=f(1)=f(-1);

If sin(c)=0, then c=n*pi, where n is an integer. This would not work because cos(n*pi) does not equal zero, and the boundary conditions would not be true.

> assume(sin(c)=NumeralNonZero);

> B:=solve(eqn5,B);

> Theta(eta,tau);

> eqn6:=0=f(1);

> eqn7:=0=f(-1);

Let's see all of the solutions.

> _EnvAllSolutions:='true';

> solve({eqn6,eqn7},{c,C});

This is the eigenvalue. Cosine of this eigenvalue is the eigenfunction.

> c:=Pi*(n+1/2); assume(n,integer);

The eigenfunctions form a basis for the general solution. Thus we may multiply them together.

> simplify(Theta(eta,tau));

A and C may actually be different for each n. The most general solution of this form is obtained by adding the solutions of the form of Eq. 11.1-25 for all integral n from n= -infinity to n=+infinity. However, we will define a integration constant that allows us to sum from n=0 to n=+infinity. We will call it Dn as in the book (see Eq 11.1-26).

> Theta:=(eta,tau)->Dn*cos(1/2*Pi*(2*n+1)*eta)*exp(-1/4*Pi^2*(2*n+1)^2*tau);

This is Eq. 11.1-26

> Solution:=(eta,tau,N)->Sum(Theta(eta,tau),n=0..N);

This is Eq. 11.1-27. We get this by using the initial condition that tau=0 when Theta=1.

> eqn:=1=Solution(eta,0,N);

This is Eq. 11.1-29

> eqn2:=int(lhs(eqn)*cos((m+1/2)*Pi*eta),eta=-1..1)=int(rhs(eqn)*cos(((m+1/2)*Pi*eta)),eta=-1..1);

We will evaluate the right hand side by taking the summation outside of the integral.

> Right:=simplify(Sum(int(cos((m+1/2)*Pi*eta)*Dn*cos(1/2*Pi*(2*n+1)*eta),eta=-1..1),n=0..N));

When the integration is done, all the integrals in the summation on the right are identically zero, as we can see above, except when n=m. This is because all of the eigenfunctions are orthogonal. When we take the inner product of two orthogonal functions, it is zero. The only time that you do not have two orthogonal functions is when n=m. Thus we can drop the summation and take the limit of the rest as n approaches m. We drop the Dn term also so that we can solve for it later.

> limit((sin(Pi*(m-n))*m+sin(Pi*(m-n))+sin(Pi*(m-n))*n+sin(Pi*(m+1+n))*m-sin(Pi*(m+1+n))*n)/(m^2+m-n-n^2)/Pi,n=m);

This limit represents the right hand side of Eq. 11.1-29. To solve for Dm, we divide the left hand side of eqn2 by this limit.

> Dm:=simplify(lhs(eqn2)/%);

We now substitute a trig identity for the numerator and take advantage of the fact that sin(Pi*m) is zero and cos(Pi*m) is (-1)^m. We want Dn back again for our final solution. This is exactly the same as Eq. 11.1-30 by inspection.

> Dn:=subs({m=n,sin(Pi*m)=0,cos(Pi*m)=(-1)^n,sin(1/2*(2*m+1)*Pi)=sin(m*Pi)*cos(Pi/2)+cos(m*Pi)*sin(Pi/2)},Dm);

This is Eq. 11.1-31 and matches exactly.

> Solution(eta,tau,N);

> plot({Solution(eta,1,100),Solution(eta,.4,100),Solution(eta,.04,100)},eta=0..1);

These are a couple of plots of the temperature profile at tau=1, 0.4, and 0.04. They match with the book's plots.