• Packages
• Introduction
• Definitions
• Dirichlet Boundary Conditions
  • f(x) = sin(Pi*x/L)+2*sin(4*Pi*x/L);
  • f(x) = 1;
• Newmann Boundary Conditions
  • f(x) = -2*sin(Pi*x/L);
• References
• Disclaimer

Partial Differential Equations

The Heat Equation: Separation of variables and Fourier series

Anton Dzhamay

Department of Mathematics

The University of Michigan

Ann Arbor, MI 48109

wPage: http://www.math.lsa.umich.edu/~adzham

email: adzham@umich.edu

Copyright  2004  by Anton Dzhamay

All rights reserved

Packages

Some packages that we use in this worksheet:

>	restart: with(plottools): with(plots):

Warning, the names arrow and changecoords have been redefined

Introduction

In this worksheet we consider the one-dimensional heat equation   describint the evolution of temperature    inside the homogeneous metal rod. We consider examples with homogeneous Dirichlet (   ,   ) and Newmann (   ,   )  boundary conditions and various  

initial profiles  . Say, we want to solve the problem with homogeneous Dirichlet boundary conditions. Using separation of variables we can get an infinite family of particular solutions of the form     . Then we take a linear combination    of such solutions with the coefficients chosen in such a way that at    we get the initial profile   .   

Definitions

>	assume(n,integer); assume(m,integer);

Shorthand notation for basic functions

>	s:=(x,n)->sin(n*Pi*x/L):s(x,n); c:=(x,n)->cos(n*Pi*x/L):c(x,n);

Particular solutions

>	uSp:=(x,t,n)->s(x,n)*exp(-(n*Pi/L)^2*k*t):uSp(x,t,n); uCp:=(x,t,n)->c(x,n)*exp(-(n*Pi/L)^2*k*t):uCp(x,t,n);

Fourier sine coefficients for    on the interval  

>	B:=proc(expr,var,n)         simplify(int(expr*s(var,n),var=0..L)/int(s(var,n)*s(var,n),var=0..L)); end proc:B(f(x),x,n);

Fourier cosine coefficients for    on the interval    (note that the formulas are different for    and   )

>	A:=proc(expr,var,n)         simplify(int(expr*c(var,n),var=0..L)/int(c(var,n)*c(var,n),var=0..L)); end proc:A(f(x),x,0);A(f(x),x,n);

Fourier sine series and Fourier sine polynomial  for  [The subtle difference here is that sometimes series (that uses sum ) has troubles with division by zero. The polynomial (that uses add ) does not have this problem, but on the other hand can not evaluate symbolic sums]

>	FPs:=proc(expr,var,n)         add(B(expr,var,m)*s(var,m),m=1..n); end proc: FSs:=proc(expr,var,n)         sum(B(expr,var,m)*s(var,m),m=1..n); end proc:FSs(f(x),x,infinity);

Fourier cosine series and Fourier cosine polynomial for  

>	FPc:=proc(expr,var,n)         A(expr,var,0)+add(A(expr,var,m)*c(var,m),m=1..n); end proc: FSc:=proc(expr,var,n)         A(expr,var,0)+sum(A(expr,var,m)*c(var,m),m=1..n); end proc:FSc(f(x),x,infinity);

Fourier polynomial (and series) solution for homogeneous Dirichlet boundary conditions

>	uPs:=proc(expr,xvar,tvar,n)         add(B(expr,xvar,m)*s(xvar,m)*exp(-(Pi*m/L)^2*k*tvar),m=1..n); end proc: uSs:=proc(expr,xvar,tvar,n)         sum(B(expr,xvar,m)*s(xvar,m)*exp(-(Pi*m/L)^2*k*tvar),m=1..n); end proc:uSs(f(x),x,t,n);

Fourier polynomial (and series) solution for homogeneous Neumann boundary conditions

>	uPc:=proc(expr,xvar,tvar,n)         A(expr,xvar,0)+add(A(expr,xvar,m)*c(xvar,m)*exp(-(Pi*m/L)^2*k*tvar),m=1..n); end proc: uSc:=proc(expr,xvar,tvar,n)         A(expr,xvar,0)+sum(A(expr,xvar,m)*c(xvar,m)*exp(-(Pi*m/L)^2*k*tvar),m=1..n); end proc:uSc(f(x),x,t,infinity);

>

Dirichlet Boundary Conditions

 

In this example we take    and   .

>	L:=1: k:=1:

Since  is just the sum of two basic functions for  and , so by the principle of superposition  is the required solution.

>	u:=(x,t)->uSp(x,t,1)+2*uSp(x,t,4):u(x,t);

This plot shows our solution and the corresponding particular solutions. Looking at their evolution in time we see that due to faster exponential decay the particular solution    goes to zero much faster than    and so when    is large enough the solution is well-approximated by just the first particular solution.

>	display([   seq(     plot({u(x,0.005*t),uSp(x,0.005*t,1),2*uSp(x,0.005*t,4)},x=0..L,      color=[red,blue,green],thickness=2)                   ,t=0..50) ],insequence=true);

And this is how the graph of our solution looks like from the top (we extend it a bit in the   -direction and also scale the   -values by the coefficient of 1/10 to get a better picture). Points are colored according to the temperature, and we can see how  quickly the initial sharp variations in the temperature (which correspond to large   ) disappear:

>	display([seq(plot3d(u(x,0.005*t)/10, x=0..L,y=0..L/8,shading=zhue,style=patchnogrid),t=0..50)],insequence=true, scaling=constrained,axes=boxed,orientation=[-90,0],tickmarks=[1,1,1]);

Now let's restore    and    back to symbols

>	L:='L': k:='k':

 

This time we'll need the full power of Fourier series.

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

The Fourier sine coefficients of    are given by

>	B(f(x),x,n);

Here are the first few of them evaluated explicitly (note that all odd coefficients are zero in this case).

>	seq(B[n]=B(f(x),x,n),n=1..10);

So the Fourier sine series for    is

>	FSs(f(x),x,infinity);

And this is an approximation whith four first non-zero terms:

>	FSs(f(x),x,8);

For graphing we again choose some values for    and   .

>	L:=1: k:=1:

Sketching Fourier sine series approximations we see that away from the boundary points Fourier approximation is getting closer and closer to our function

>	display([ plot(1,x=0..L,color=red,thickness=2), display([seq(plot(FSs(f(x),x,2*r+1),x=0..L,color=blue,thickness=2, title=sprintf("n=%d Fourier Approximation",2*r+1)),r=0..20)],insequence=true) ]);

Sketching the Fourier approximation with    helps us see it even better. Near    we see that Fourier approximation goes to about 1.18. This is known as the Gibbs phenomena for the Fourier series.

>	plot([f(x),FSs(f(x),x,300)],x=0..L,color=[red,blue], thickness=[3,1],numpoints=1000,scaling=constrained);

Note, however, that outside of the interval    the situation is quite different. Fourier series of    converges not to    but to its odd periodic extension .

>	plot([f(x),FSs(f(x),x,100)],x=-3*L..3*L,color=[red,blue], thickness=[2,1],numpoints=1000,scaling=constrained);

Looking at the   -dynamics of the Fourier polynomial of order 15 (red below) and its different term we again see that higher harmonics (green) have amplityde decaying faster than the first harmonic (principal mode, blue):

>	u:=(x,t)->uPs(f(x),x,t,15): display([   seq(     plot([u(x,0.005*t),seq(B(f(x),x,2*r+1)*uSp(x,0.005*t,2*r+1),r=0..7)],x=0..L,     color=[red,blue,seq(green,r=1..7)],thickness=2)           ,t=0..40) ],insequence=true);

Another way to see it is to look at the "energies" of different mode. We define the energy of    to be its   -norm,   :

>	eMode:=(n,t)->int(uSp(x,t,n)^2,x=0..L):eMode(n,t);

 

This is the relative decay speed of different modes:

>	display([seq(  display([seq(disk([0.2*r,eMode(2*r+1,0.005*t)],0.02,color=red),r=0..7)]) ,t=0..40)],insequence=true,scaling=constrained,view=[-0.1..1.5,0..1]);

And this is their actual contribution to the initial profile - note that the first mode makes the "bulk" of the initial energy, and also decays the slowest!

>	display([seq(  display([seq(disk([0.2*r,B(f(x),x,2*r+1)^2*eMode(2*r+1,0.005*t)],0.02,color=red),r=0..7)]) ,t=0..40)],insequence=true,scaling=constrained,view=[-0.1..1.5,0..1]);

Clean-up:

>	L:='L': k:='k': f:='f':

Newmann Boundary Conditions

 

>	f:=x->-2*sin(Pi*x/L):f(x);

Fourier cosine coefficients:

>	A(f(x),x,n);

Explicitly:

>	seq(A[n]=A(f(x),x,n),n=0..10);

This is how the Fourier polynomial of order 4 looks like:

>	FPc(f(x),x,4);

For graphing we again take , .

>	L:=1: k:=1:

This function plots the graph of  and its -th approximation on the interval .

>	approxplot:=k->plot([f(x),FPc(f(x),x,k)],x=0..L, color=[red,blue],title=sprintf("%d-th approximation",k),view=-2.3..0.3, legend=["Initial profile",sprintf("%d-th approximation",k)], linestyle=[SOLID,DASHDOT]):

Here is a sequence of such approximation. Notice the convergence:

>	display([seq(approxplot(5*i),i=0..6)],insequence=true);

Note, however, that outside of the  interval the function  and its Fourier polynomial are quite different:

>	plot([f(x),FPc(f(x),x,15)],x=-3*L..3*L,color=[red,blue]);

The function below sketches the initial profile and the Fourier polynomial approximation to the solution of order  at time .

>	timeplot:=(t,k)->plot([f(x),uPc(f(x),x,t,k)],x=0..L, color=[red,blue],title=sprintf("%d-th FS approximation at time t=%2.2f",k,t),view=-2.3..0.3,legend=["Initial profile", sprintf("Profile at time t=%2.2f",t)],linestyle=[SOLID,DASHDOT]):

We can see that the solution converges to the average of the initial profile .

>	display([timeplot(0,20),timeplot(0.01,20),timeplot(0.05,20),timeplot(1,20)],insequence=true);

>

References

• Walter A. Strauss, Partial Differential Equations: An Introduction, Wiley, 1992
• Richard Haberman, Elementary Applied Partial Differential Equations, 3rd edition, Prentice Hall

Disclaimer

"While every effort has been made to validate the solutions in this worksheet, Waterloo Maple Inc. and the contributors are not responsible for any errors contained and are not liable for any damages resulting from the use of this material."