#############################################################

# The Heat Equation

# Jim Carlson 
# http://www.math.utah.edu/~carlson)
# December 14, 1995

# The code below solves a discrete version of the one-dimensional
# heat equation.  The initial temperature state is stored in 
# the array n-element array u.  One can use zero(u) to set
# this array to zero, and then set individual elmenents by
# statements such as u[5] := 1;  The procedure call
# update(u, uu) computes the next temperature state and stores
# it in the array uu.  See the examples below for 
# more information, including directions for making movies
# using the film and show commands.

# The main idea: the temperature in a given cell at time
# t+1  tends to be closer to the nearby average temperature
# than it was at time  t.  Denote that temperatue in cell
# k  at time  t  by  u[k], and the temperature at time  t+1
# in that cell by  uu[k].  One measure of the local average
# is   
#            (u[k-1] + u[k+1])/2
#
# One way of expressing the tendency to move towards the 
# local average is to say
#
#     uu[k+1] := u[k] + c( (u[k-1] + u[k+1])/2 - u[k] )
#
# where c is a positive constant of proportionality.
# We have used c = 1/2 in the code below.

# To add:

# References: .... to be found ....

# Note: (on numerical instability)

# Note: (on differential and difference equations)


#############################################################

# Set up variables:

n := 11;
u := array(1..n):
uu := array(1..n):

#############################################################

# Define procedures:

dim := a -> op(2,op(2,eval(a))):

zero := proc(s)
   local i;
   for i from 1 to dim(s)  do
     s[i] := 0.0;
   od;
   s;
end:

update := proc( s, ss )
   local i, n;
   n := dim(s);
   ss[1] := s[1]; ss[n] := s[n];
   for i from 2 to n-1 do
     ss[i] := s[i]/2.0 + (s[i-1] + s[i+1])/4.0;
   od;
   ss;
end:

cp := proc( target, source )
    local i;
    for i from 1 to dim(source) do
       target[i] := source[i];
     od;
end:

plotstate := proc( s )
   local i, g;
   g := array[1..dim(s)];
   for i from 1 to dim(s) do
      g[i] := [ i-1, s[i] ];
   od;
   # print(t);
   plot(convert(g,list));
end:

film := proc( u, T )
    local i, gu;
    gu := array(0..T);
    gu[0] := plotstate(u);
	for i from 1 to T do
	   update( u, uu );
	   gu[i] := plotstate(uu);
	   cp(u,uu);
	od:
    gu;
end:

show := gu -> display( convert(gu,list), insequence=true);


#############################################################

# Example 1. (Dirac delta function)

zero(u);
u[6] := 1;

print(u);
plotstate(u);

update(u,uu);

print(u);
plotstate(u);

#############################################################

# Example 2.  Now we will use a loop.
# cp(u,uu) copies the contents of uu into u
# T is the number of states computed after
# the initial state.

zero(u);
u[6] := 1;
T := 4;

Digits := 3;

for i from 1 to T do
   update( u, uu );
   print(uu);
   cp(u,uu);
od:

Digits := 10;

#############################################################

# Example 3. A film.
# Be sure to load the plots package.


zero( u );
u[6] := 1;
T := 10;

with(plots):


show( film(u,T) );

#############################################################