The heat equation and Fourier series

We will study the solution to the heat equation of a wire of length L. The body of the wire is insulated but the ends may or may not be (depending on the boundary conditions we shall set). In the particular example below, the ends are not be insulated but instead will be kept at 0 degrees at all times.

sm3110_heat1.mws,wdj,2-4-98

> with(plots):

> heat:=a->a*diff(u(x,t),x$2)=diff(u(x,t),t):
heat(a);
This is the partial diffrential equation governing the temperature function u(x,t). We wil plot the exact expression for u(x,t) (due to Fourier and which uses Fourier series) below. The constant a>0 depends on the conductivity of the wire.

We shall take L=Pi and a=1 in the example below.

The program below is gives the n-th sine series coefficient.

> sine_series_coeff:=proc(f,L,n::integer)
#f is a function of x, 0<x<L
#n is the index of b_n
#L is half the period
local x,b;
b:=2*int(f(x)*sin(n*Pi*x/L),x=0..L)/L;
RETURN(b);
end;

We first specify an initial temperature for a wire.

> f:=x->piecewise(x<Pi/2,1,x>Pi/2,2);

> plot(f(x),x=0..Pi,y=0..2);

Some partial sine series of this function are plotted below to illustrate how they approximate the function. The reason for using sine series to solve a zero ends heat equation is because the heat equation is much easier to solve when the initial temperature is a sine or a cosine function and the sine function matches the zero ends boundary conditions.

> S2:=add(sine_series_coeff(f,Pi,n)*sin(n*x),n=1..2);
S5:=add(sine_series_coeff(f,Pi,n)*sin(n*x),n=1..5):
S10:=add(sine_series_coeff(f,Pi,n)*sin(n*x),n=1..10):
S15:=add(sine_series_coeff(f,Pi,n)*sin(n*x),n=1..15):
plot([f(x),S2,S5,S10,S15],x=0..Pi);

The program below will return the sum of the first m terms of the solution to a zero ends heat equation using the Fourier series method.

> heat:=proc(f,L,a,m::integer,x,t)
local n,s;
s:=add((sine_series_coeff(f,L,n)*sin(n*Pi*x/L)*exp(-a*(n*Pi/L)^2*t)),n=1..m):
RETURN(s);
end;

>

> heat(f,Pi,1,2,x,t);
The first 2 terms of the solution. This approximates u(x,t) fairly well if t is large (say, t>2).

> P1:=plot(heat(f,Pi,1,30,x,0),x=0..Pi,y=0..2.2):
P2:=plot(heat(f,Pi,1,30,x,.05),x=0..Pi,y=0..2.2):
P3:=plot(heat(f,Pi,1,30,x,.1),x=0..Pi,y=0..2.2):
P4:=plot(heat(f,Pi,1,30,x,.2),x=0..Pi,y=0..2.2):
P5:=plot(heat(f,Pi,1,30,x,1),x=0..Pi,y=0..2.2):
P6:=plot(heat(f,Pi,1,30,x,2),x=0..Pi,y=0..2.2):
P7:=plot(heat(f,Pi,1,30,x,4),x=0..Pi,y=0..2.2):
P8:=plot(heat(f,Pi,1,30,x,6),x=0..Pi,y=0..2.2):
P9:=plot(heat(f,Pi,1,30,x,8),x=0..Pi,y=0..2.2):
P10:=plot(evalf(heat(f,Pi,1,30,x,10)),x=0..Pi,y=0..2.2):

> display([seq(P.i,i=1..10)]);

> animate(heat(f,Pi,1,50,x,t),x=0..Pi,t=0..10,title=`heat equation`);