29hotRect.html

Partial Differential Equations PowerTool

by Dr. Jim Herod

Section 7.1: The Heat Equation on a Rectangle

Maple Packages for Section 7.1

> restart;

>

We consider the diffusion of heat into a long beam with cross section a rectangle. The supposition that the beam is "long" is to produce the mathematical idea that heat diffusing to a point comes essentially only from the sides and that the ends are so far away that heat coming from above or below can be ignored. This is the same problem as considering heat diffusion in a thin rectangular plate that is insulated at the top and bottom. These two problems lead to an equation such as

with boundary conditions

u( t, x, 0) = , u( t, x, b) =

u( t, 0, y) = , u( t, a, 0) =

and initial conditions

u( 0, x, y) = h( x, y).

The problem is broken into two parts: the steady state v(x,y) and the transient, w( t, x, y).

The steady-state problem is

0 =

with boundary conditions

v( t, x, 0) = , v( t, x, b) =

v( t, 0, y) = , v( t, a, 0) = .

We have discussed how to break this problem into two parts previously. They would be

0 = ,

with boundary conditions

v1( x, 0) = , v1( x, b) =

v1( 0, y) = , v1( a, 0) =

and

0 = ,

with boundary conditions

v2( x, 0) = 0, v2( x, b) = 0,

v2( 0, y) = , v2( a, 0) = .

The solution for these two problems are added to get the steady state solution

v(x, y) = v1(x, y) + v2(x, y).

We will not detail how to make this computation for the steady state solution, since we have covered how to compute solutions for Laplace's Equation on a rectangle earlier.

The transient problem is

with boundary conditions

w( t, x, 0) = 0, w( t, x, b) = 0

w( t, 0, y) = 0, w( t, a, 0) = 0

with initial condition

w( 0, x, y) = h( x, y) - v( x, y).

Think about why u(t, x, y) = w(t, x, y) + v(x, y) is the solution for the original problem.

We break this last PDE into three ODE's (separation of three variables, assuming w = X(x)Y(y)T(t)), two of which have boundary conditions.

X '' = - X, Y '' = - Y T ' = - ( ) T

X(0) = X(a) = 0 Y(0) = Y(b) = 0.

The result is a general solution of the transient equation in the form

w( t, x, y) =

The coefficients are obtained by Fourier series as

Check

We check what was purported to be a solution for the PDE.

> u:=(t,x,y)->sum(sum(A[n,m]*sin(n*Pi*x/a)*sin(m*Pi*y/b)*
exp(-(n^2/a^2+m^2/b^2)*Pi^2*t),n=1..3),m=1..3);

Here are the boundary conditions.

> u(t,0,y); u(t,a,y);
u(t,x,0); u(t,x,b);

This would be the initial condition.

> u(0,x,y);

And, here is the PDE.

> simplify(diff(u(t,x,y),t)-diff(u(t,x,y),x,x)-diff(u(t,x,y),y,y));

Example

We illustrate how an animation of such a result behaves.

Take boundary conditions as follows:

u( t, x, 0) = sin( x), u( t, x, 1) = 0

u( t, 0, y) = 0, u( t, 1, y) = 0

Take the initial condition as

u( 0, x, y) = sin( x).

By now, maybe we can do the steady state without the computation of any integrals. The formula for the steady state, and a check follows.

> restart;

> v:=(x,y)->(sinh(Pi*(1-y)))*sin(Pi*x)/sinh(Pi);

> diff(v(x,y),y,y)+diff(v(x,y),x,x);

> v(0,y);v(1,y);v(x,0);v(x,1);

The function w will be the transient solution. The initial condition for w is k, given below.

> k:=(x,y)->sin(Pi*x)-v(x,y);
k(0,y);k(1,y);k(x,0);k(x,1);

We now compute the coefficients as directed above. I set this up to do 100 double integrals. You could do more, or less depending on the speed of your computer. The point is that k(x,y) is not continuous on the closed rectangle making the domain of u. Thus, we will not have a good approximation for u(0, x, y).

>    for n from 1 to 10 do
  for m from 1 to 10 do int(int(k(x,y)*sin(n*Pi*x)*sin(m*Pi*y),x=0..1),y=0..1)/
     (int(sin(n*Pi*x)^2,x=0..1)*int(sin(m*Pi*y)^2,y=0..1)):
     A[n,m]:=evalf(%):
  end do:
end do:

> n:='n': m:='m':

We make the transient solution w.

> sum(sum(A[n,m]*sin(n*Pi*x)*sin(m*Pi*y)*exp(-(n^2+m^2)*Pi^2*t),
n=1..10),m=1..10):

> w:=unapply(%,(t,x,y)):

From w and v, the solution for the original problem is made.

> u:=(t,x,y)->w(t,x,y)+v(x,y);

As a check to see how accurate all this is, u( 0, x, y) should be sin( x) in this problem.

> plot3d(u(1/100,x,y),x=0..1,y=0..1,axes=NORMAL);

Finally, here is an animation of the temperature decaying to the steady state.

> plots[animate3d](u(t,x,y),x=0..1,y=0..1,t=0..1/10, frames=50);

>

Example: Speed of Diffusion.

If we want to incorporate speed of diffusion, we should change the model by including c.:

This time, we choose conditions

u( t, x, 0) = 100, u( t, x, 1) = 100

u( t, 0, y) = 100, u( t, 1, y) = 100

and initial conditions

u( 0, x, y) = 0.

We ask, for various values of c, what is the value of t such that u(t, 1/2, 1/2) = 50.

Is it clear that the steady state solutions is 100? And, that general solution is

u( t, x, y) = + 100,

where

= -400 .

>    for n from 1 to 20 do
  for m from 1 to 20 do
     A[n,m]:=-400*(cos(n*Pi)-1)/(n*Pi)*(cos(m*Pi)-1)/(m*Pi):
  end do:
end do:
n:='n': m:='m':

> sum(sum(A[n,m]*sin(n*Pi*x)*sin(m*Pi*y)*exp(-c*(n^2+m^2)*Pi^2*t),
n=1..20),m=1..20):

We write this u as a function of t, x, y, and also c.

> u:=unapply(%,(c,t,x,y))+100:

Here, we plot the initial value for c = 1. The purpose is to see how well we have fit the initial condition.

> plot3d(u(1,0,x,y),x=0..1,y=0..1,axes=NORMAL, numpoints=1000);

>    simplify(diff(u(c,t,x,y),t)-
         c*diff(u(c,t,x,y),x,x)-
               c*diff(u(c,t,x,y),y,y));

Here, we check the boundary conditions.

> u(1,t,x,0);u(1,t,x,1);u(1,t,0,y);u(1,t,1,y);

This is an animation so that you can see how the solution has limit the steady state.

> plots[animate3d](u(1,t,x,y),x=0..1,y=0..1,t=0..1/10,axes=NORMAL, frames=50);

>

The goal now is to choose a sequence of c 's and, for each c, to find t so that the point [1/2, 1/2] is at temperature half way between 0 and 100. We restrict solutions for t to being in the interval [0, 5].

> for p from 1 to 10 do
s[p]:=fsolve(u(p/10,t,1/2,1/2)=50.,t,0..5);
od;
p:='p';

Here, we plot the values of t found above. The purpose is to see the kind of relationship there is between the speed of diffusion and c.

> plots[pointplot]([seq([p,s[p]],p=1..10)]);

>

What this graph shows is that as c increases the time for the middle point to reach the temperate 50 decreases. The decrease seems like an exponential decay. Can you justify that in the mathematics, or in an exponential fit to the data?

EMAIL: herod@math.gatech.edu or jherod@tds.net

URL: http://www.math.gatech.edu/~herod