17convec.html

Partial Differential Equations PowerTool

by Dr. Jim Herod

Section 4.4: Convection Across Boundaries

Maple Packages for Section 4.4

> restart;

> with(plots):

Warning, the name changecoords has been redefined

Suppose that we take a well insulated rod -- sides and ends insulated-- which has initial heat distribution as indicated by the function

4 x (1 - x ) + 2, 0 < x < 1.

We can see that this rod will be hottest in the middle, and that as the heat diffuses, a uniform distribution of heat will be achieved with total heat the same as the total heat of the original distribution. After all, no heat has been lost or gained in this well-insulated rod. Here is the definition of the initial distribution of heat, and the value of the total heat in the system.

> f:=x->4*x*(1-x)+2;

> int(f(x),x=0..1);

>

How to model the redistribution of heat was discussed in Section 4.3. Consider this alternate problem. Keep the rod well insulated, except at one end. There heat escapes or gathers according as to whether it is more or less than the outside temperature equal to 2. More specifically, We suppose

The PDE: ,

The Boundary Conditions: (t, 0) = 0 and (t,1) = - ( u (t, 1) - 2).

The Initial Condition: u (0, x) = 4 x (1-x) + 2.

Now this problem is more interesting. From the equation for the distribution for the initial temperature, we can calculate that the temperature is 2 at the two ends and 3 in the exact middle. The heat begins to diffuse, causing temperatures to rise or fall along the rod. As soon as the temperature goes above 2 on the right, the rod begins to cool by heat loss at that end. Heat will continue to leak out until the temperature over the entire rod is 2. But what happens at the left side. Surely heat will start to spread there from the middle. We calculated in the lines above that with no leakage out the right side or anywhere else, the temperature at the left would rise to 8/3. How high will it go now with this leakage out the right side? Will it go above 2, and then settle back down? And what about the hot spot? Will it move left as heat leaks out the right side, or will it drop straight down?

Maybe you agree this problem will be interesting to consider.

The PDE is homogeneous, but the boundary conditions are not. Homogeneous boundary conditions would be

The Boundary Conditions: (t, 0) = 0 and (t,1) = - w(t, 1).

Thus, we need to find a particular solution, or what we call the steady state solution, for the original problem. Is it clear that the steady state solution is

v(x) = 2?

That solution satisfies the PDE and the non-homogeneous boundary conditions. The u which satisfied the original problem is related to the transient solution w

by u = w + 2.

Now, we try to find the general solution for the PDE with homogeneous boundary conditions.

The general solution for the homogeneous problem

By now, we know to perform what was called Step 1 in Section 4.1: using separation of variables create two ordinary differential equations with boundary conditions whose solutions will determine the solution for the partial differential equation:

X '' = X, with X'(0) = 0 and X '(1) +X(1) = 0,

and

T ' = T.

Step 2 was to solve the two ordinary differential equations with boundary conditions. The constant must be negative, call it . Solutions for the differential equation

X '' = X

are A sin( x) + B cos( x).

We invoke the boundary conditions.

> X:=x->A*sin(lambda*x)+B*cos(lambda*x);

The left boundary condition is used first.

> D(X)(0)=0;

This condition implies A = 0.

> A:=0;

We now use the remaining boundary condition.

> D(X)(1)+X(1)=0;

This implies that tan( ) = .

We must find 's that satisfy this requirement. Then, = and the eigenfunction corresponding to this eigenvalue is cos( x). We have only to find the 's.

Is it now clear that there is going to be no nice, closed form solution for these 's? We must find them numerically. To help in doing this, we look at where the graph of tan(x) crosses the graph of 1/x.

> plot([tan(x),1/x],x=0..20,y=0..3/2,discont=true,color=[RED,BLACK]);

We see there are roots between all the places where tan(x) is zero and where it is undefined.

Maybe seven will be enough.

> for n from 1 to 7 do
lambda[n]:=fsolve(tan(x)=1/x,x,(n-1)*Pi..(2*n-1)*Pi/2);
od;

>

We now have approximations for seven of the eigenvalues. We are ready to take Step 3: construct the general solution for the partial differential equation. We have

w(t, x) = .

We go back to solve the original problem. We need the general solution plus the particular solutions. These add together to give the solution to this problem.

The solution for the original problem

Now, we find the c's by using the initial condition: f(x) = w(0, x) + 2

or f ( x ) -2 = .

We compute coefficients.

>    for n from 1 to 7 do
   c[n]:=int((f(x)-2)*cos(lambda[n]*x), x=0..1)/
                    int(cos(lambda[n]*x)^2,x=0..1);
od;
n:='n':

We check to see how close this is to our f .

> plot([f(x),2+sum(c[n]*cos(lambda[n]*x),n=1..7)],x=0..1);

We now define the solution, check to see that it really is a solution, and plot its graph.

> u:=(t,x)->2+sum(c[n]*exp(-lambda[n]^2*t)*cos(lambda[n]*x),n=1..7);

Check the solution

We check the differential equation and boundary conditions.

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

> D[2](u)(t,0); D[2](u)(t,1)+(u(t,1)-2);

The above answers might be a surprise. They are not all zeros. Not to worry.
The 's we computed above were only an approximation for the real ones. (Remember that we used fsolve?) So, the numbers computed just above that we thought would be zero are approximations for zero. Of course! Just look at their magnitude.

Finally, we observe how good a fit we have for the initial value.

> plot([f(x),u(0,x)],x=0..1);

>

Graph the solution

> plot3d(u(t,x),x=0..1,t=0..3,axes=NORMAL,orientation=[-20,60]);

Does this graph have the properties we predicted? Does it appear that the graph changes from being the same shape at t = 0 as the graph of f ( x ) and that the graph changes toward being constant 2 as t evolves.

We are satisfied that we have an approximation for the solution for the original problem.

In order to watch the maximum move from one side to the other, we animate u ( t , x ) as t increases. Using the cursor, touch the graph of the animation which comes next. Use the Move To The Next Frame button to watch how the maximum temperature moves to the left.

> animate(u(t,x),x=0..1,t=0..1,frames=100,view=[0..1,1.5..3.5]);

We ask: At what time would the left end point and right end point have maximum temperatures? This graph gives some answer. The red graph is the left endpoint temperature and the green is the right endpoint temperature.

> plot([u(t,0),u(t,1)],t=0..2,y=0..3,color=[red,green]);

>

This problem was important for it led to an eigenvalue problem for which we could not find the eigenvalues in closed form. It was necessary to use numerical methods to generate these eigenvalues. From an applications perspective, the boundary conditions are more realistic for it is unlikely that the end points can be held fixed or perfectly insulated.

Unassisted Maple

We have seen how Maple suggests separation of variables in previous sections for the diffusion equation and variations. In this section, we allow Maple freedom on a different problem.

Perhaps here the user felt uncomfortable with having to make approximations for where solutions for the equation

tan( ) = 1/

should lie. After all, if we wanted 20 or 200 such 's, we would not want to go searching for the intervals in which to seek solutions. The question we address is how to let Maple find these without human assistance.

The answer is to ask the right question. This question is suggested by the following graphs of 1/x and tan(x).

> plot([tan(x),1/x],x=0..20,y=-2..2, discont=true);

This graph makes clear the obvious. There is a solution for the equation tan(x) = 1/x between successive odd multiples of /2. Thus, we let Maple determine these values unassisted by our searching out specific intervals for each solution.

> for n from 0 to 7 do
mu[n]:=fsolve(tan(x)=1/x,x,max(0,(2*n-1)*Pi/2)..(2*n+1)*Pi/2);

> od;

>

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

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