30insul.html

Partial Differential Equations PowerTool

by Dr. Jim Herod

Section 7.2: Two Dimensional Diffusion with Neumann Boundary Conditions

Maple Packages for Section 7.2

> restart:

> with(plots):

Warning, the name changecoords has been redefined

>

In Lecture 32, we considered boundary conditions that are called Dirichlet Boundary Conditions . Here, we will work a problem with Neumann Conditions corresponding to an insulated surface in the diffusion of heat.

We illustrate this type problem with a simple example.

Consider the heat conduction problem in a rectangle where there is insulation on all boundaries and the initial condition is as follows:

Insulated boundaries: (t,x,0) = (t,x,1)= 0, (t,0,y) = (t,1,y)= 0.

Initial condition: u(0, x, y) = 100 x (1-x) y (1-y).

Here are the tasks we will accomplish:

(A) Give three ODE's, together with boundary conditions, whose solutions lead to a solution for this problem.

(B) Give the solution for the three ODE's.

(D) Approximate the particular solution that goes with this initial condition.

(E) Animate a graph of the solution.

(F) Compute the total heat for the system for all time.

(A.) Here are three ODE's and boundary conditions:

X '' + X = 0, X '(0) = X '(1) = 0.

Y '' + Y = 0, Y '(0) = Y '(1) = 0.

T ' + ( ) T = 0.

(B.) Solutions for these equations are

X(x) = cos( n x), Y(y) = cos(m y), T(t) = exp( - ( ) t )

(C.) The general solution for the pde is

u(t,x,y) = .

This is so important I will check an approximation of the sum.

Here defines u.

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

Here checks the boundary conditions.

> D[3](u)(t,x,0);
D[3](u)(t,x,1);
D[2](u)(t,0,y);
D[2](u)(t,1,y);

Here checks the PDE.

> D[1](u)(t,x,y)-D[2,2](u)(t,x,y)-D[3,3](u)(t,x,y);

(D.) We compute several terms of the particular solution.

> 100*x*(1-x)*y*(1-y)=u(0,x,y):

>    for n from 0 to 10 do
  for m from 0 to 10 do
    A[n,m]:=100*int(x*(1-x)*cos(n*Pi*x),x=0..1) /
            int(cos(n*Pi*x)^2,x=0..1)
            * int(y*(1-y) *cos(m*Pi*y),y=0..1)/int(cos(m*Pi*y)^2,y=0..1);
  end do:
end do:
n:='n': m:='m':

>

In order to check my solution, I compare the graph of the initial value with the initial function.

> plot3d(u(0,x,y),x=0..1,y=0..1,axes=normal, shading=zhue);

> plot3d(100*x*(1-x)*y*(1-y),x=0..1,y=0..1,axes=normal, shading=zhue);

Note that there is a small error along the boundaries. This error could be made smaller by computing more than the 121 coefficients that we have computed above.

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

(E.) Here is an animation of the solution.

> animate3d(u(t,x,y),x=0..1,y=0..1,t=0..1/8,axes=normal, frames=30);

(F.) We compute the total heat and verify that it is the same as the initial total heat.

> int(int(u(t,x,y),x=0..1),y=0..1);
TotHeat:=int(int(100*x*(1-x)*y*(1-y),x=0..1),y=0..1);

>

Finally, we recall how to make changes in case the boundary conditions were not homogeneous. A quick review for working such a problem seems appropriate.

Non homogeneous Neumann Type Boundary Conditions.

Suppose the partial differential equation and boundary conditions have the form

Boundary Conditions: (t,x,0) = 0, (t,x,1)= -1, (t,0,y) = 0, (t,1,y)= 1.

Initial condition: u(0, x, y) = x (1-x) y (1-y) + ( ) / 2.

Note changes from the previous problem:

The boundary conditions have been changed to model heat coming in at the top and going out the right side. The initial conditions have been changed so that the results computed at the beginning of this worksheet can be used here.

The pattern for working problems of this type have been set. First, we solve the steady state problem. The steady state problem will have this form:

PDE: 0 = v,

Boundary Conditions: (t,x,0) = 0, (t,x,1)= -1, (t,0,y) = 0, (t,1,y)= 1.

We have addressed this problem previously. We found the following solution which we define and check here.

> v:=(x,y)->(x^2-y^2)/2;

> D[2](v)(x,0);
D[2](v)(x,1);
D[1](v)(0,y);
D[1](v)(1,y);

> D[1,1](v)(t,x,y)+D[2,2](v)(t,x,y);

>

Next we find the transient solution which we call w . To do this we use homogeneous boundary conditions and an initial condition to reflect that the solution u for the original problem will be defined by

u (t, x, y) = w (t, x, y) + v (x, y).

The initial conditions of this non-homogeneous problem have been created so that the transient solution is exactly the problem we worked above: w (t, x, y) = u (t, x, y) - v(x,y) . We use these results here.

> w:=(t,x,y)->sum(sum(A[m,n]*cos(n*Pi*x)*cos(m*Pi*y)*
exp(-(n^2+m^2)*Pi^2*t),n=0..4),m=0..4);

The solution for the problem in this section should now be the sum of w and v . We check this.

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

We check the boundary conditions.

> D[3](u)(t,x,0);
D[3](u)(t,x,1);
D[2](u)(t,0,y);
D[2](u)(t,1,y);

We check the PDE.

> D[1](u)(t,x,y)-D[2,2](u)(t,x,y)-D[3,3](u)(t,x,y);

We check the initial conditions by comparing the graphs of u (0, x, y) and

100* x (1-x) y (1-y) - ( ) / 2.
In order to see that there are two graphs here, we displace the second by an amount 0.1.

> plot3d({u(0,x,y),100*x*(1-x)*y*(1-y)+(x^2-y^2)/2+0.1},
x=0..1,y=0..1,axes=normal, shading=zhue);

> plot3d({u(1,x,y),TotHeat+(x^2-y^2)/2+0.1},
x=0..1,y=0..1,axes=normal, shading=zhue);

Finally, to give us intuition, we draw the animation as the solution moves from the initial condition to the steady state condition.

> animate3d(u(t,x,y),x=0..1,y=0..1,t=0..1/20, axes=normal, frames=30, shading=zhue);

>

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

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