############################################################# # 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) ); #############################################################