Partial Differential EquationsThe Heat Equation: Separation of variables and Fourier seriesAnton DzhamayDepartment of MathematicsThe University of MichiganAnn Arbor, MI 48109wPage: http://www.math.lsa.umich.edu/~adzhamemail: adzham@umich.eduCopyright \251 2004 by Anton DzhamayAll rights reserved

Some packages that we use in this worksheet:restart: with(plottools): with(plots):

In this worksheet we consider the one-dimensional heat equation NiMvLSUlZGlmZkc2JC0lInVHNiQlInhHJSJ0R0YrKiYlImtHIiIiLUYlNiVGJ0YqRipGLg== describint the evolution of temperature NiMtJSJ1RzYkJSJ4RyUidEc= inside the homogeneous metal rod. We consider examples with homogeneous Dirichlet (NiMvLSUidUc2JCIiISUidEdGJw== , NiMvLSUidUc2JCUiTEclInRHIiIh ) and Newmann (NiMvLS0lJWRpZmZHNiQlInVHJSJ4RzYkIiIhJSJ0R0Yr , NiMvLS0lJWRpZmZHNiQlInVHJSJ4RzYkJSJMRyUidEciIiE= ) boundary conditions and various initial profiles NiMtJSJmRzYjJSJ4Rw== . 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 NiMvLSYlInVHNiMlIm5HNiQlInhHJSJ0RyomLSUkc2luRzYjKipGKCIiIiUjUGlHRjFGKkYxJSJMRyEiIkYxLSUkZXhwRzYjLCQqKCooRihGMUYyRjFGM0Y0IiIjJSJrR0YxRitGMUY0RjE= . Then we take a linear combination NiMvLSUidUc2JCUieEclInRHLSUkc3VtRzYkKiYmJSJCRzYjJSJuRyIiIi0mRiVGL0YmRjEvRjA7JSJpRyUpaW5maW5pdHlH of such solutions with the coefficients chosen in such a way that at NiMvJSJ0RyIiIQ== we get the initial profile NiMvLSUidUc2JCUieEciIiEtJSJmRzYjRic= .

assume(n,integer); assume(m,integer);Shorthand notation for basic functionss:=(x,n)->sin(n*Pi*x/L):s(x,n); c:=(x,n)->cos(n*Pi*x/L):c(x,n);Particular solutionsuSp:=(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 NiMtJSJmRzYjJSJ4Rw== on the interval NiM3JCIiISUiTEc= 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 NiMtJSJmRzYjJSJ4Rw== on the interval NiM3JCIiISUiTEc= (note that the formulas are different for NiMvJSJuRyIiIQ== and NiMyIiIhJSJuRw== )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 NiMtJSJmRzYjJSJ4Rw== [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 NiMtJSJmRzYjJSJ4Rw== 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 conditionsuPs:=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 conditionsuPc:=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);

In this example we take NiMvJSJMRyIiIg== and NiMvJSJrRyIiIg== .L:=1: k:=1:Since NiMtJSJmRzYjJSJ4Rw== is just the sum of two basic functions for NiMvJSJuRyIiIg== and NiMvJSJuRyIiJQ==, so by the principle of superposition NiMvLSUidUc2JCUieEclInRHLCYtJkYlNiMiIiJGJkYtKiYiIiNGLS0mRiU2IyIiJUYmRi1GLQ== 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 NiMtJiUidUc2IyIiJTYkJSJ4RyUidEc= goes to zero much faster than NiMtJiUidUc2IyIiIjYkJSJ4RyUidEc= and so when NiMlInRH 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 NiMlInlH -direction and also scale the NiMlInVH -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 NiMlIm5H ) 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 NiMlIkxH and NiMlImtH back to symbolsL:='L': k:='k':

This time we'll need the full power of Fourier series.f:=x->1:f(x);The Fourier sine coefficients of NiMtJSJmRzYjJSJ4Rw== are given byB(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 NiMtJSJmRzYjJSJ4Rw== isFSs(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 NiMlIkxH and NiMlImtH .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 functiondisplay([ 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 NiMvJSJuRyIkKyQ= helps us see it even better. Near NiMvJSJ4RyIiIQ== 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 NiM3JCIiISUiTEc= the situation is quite different. Fourier series of NiMtJSJmRzYjJSJ4Rw== converges not to NiMtJSJmRzYjJSJ4Rw== 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 NiMlInRH -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 NiMtJiUidUc2IyUibkc2JCUieEclInRH to be its NiMqJCUiTEciIiM= -norm, NiMtJSRpbnRHNiQqJC0mJSJ1RzYjJSJuRzYkJSJ4RyUidEciIiMvRi07IiIhJSJMRw== :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':

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 NiMvJSJMRyIiIg==, NiMvJSJrRyIiIg==.L:=1: k:=1:This function plots the graph of NiMtJSJmRzYjJSJ4Rw== and its NiMlImtH-th approximation on the interval NiM3JCIiISUiTEc=.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 NiM3JCIiISUiTEc= interval the function NiMtJSJmRzYjJSJ4Rw== 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 NiMlImtH at time NiMlInRH.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 NiMtJSJmRzYjJSJ4Rw==.display([timeplot(0,20),timeplot(0.01,20),timeplot(0.05,20),timeplot(1,20)],insequence=true);

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

"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."