A substantial project; the heat equation

In this project we will hopefully accomplish four goals;
1. Study simple algorithms for solving second order partial differential equations,
2. execute the algorithms with extensive use of array variables, pointers, and functions that process arrays,
3. begin a study of computer graphics and animation using third-party graphics libraries,
4. manage a programming project with make.

The heat equation

Heat flow is a diffusive process that obeys the differential equation

2 ∂-T-(x,-t)-= D ∂T-(x,t) ∂x2 ∂t
in the sense that the temperature at point x at time t is the function T(x,t) which is the solution to this equation, basically Newton’s law of cooling with overall energy conservation.
The equation is solved numerically by discretizing space and time
T (x,t) = T(n, m), x = n ⋅ dx, t = m ⋅ dt
to obtain a difference equation
T-(x-+-dx,t) --2T-(x,t) +-T-(x---dx,t) = D T(x,t-+-dt) --T-(x,t) (dx2 )2 dt
or
( --8dt--) (--4-dt--)( ) T (x,t + dt) = 1 - D(dx)2 T(x,t) + D(dx)2) T(x + dx, t) + T (x - dx,t)
We will store the temperature variable in an array, T[n] = T(n ⋅ dx,t) and use this formula to obtain an updated array Tnew[n] = T(n ⋅ dx,t + dt). This is very typical of numerical schemes to solve PDEs.

Data storage and update routines

We will write two routines that will a; find Tnew[n] from T[n], and b; swap Tnew[n] ↔ T[n] to prepare for the next iteration. These routines will be void and will accept pointers to the arrays. We will offset the arrays so that space is divided into 2N + 1 discrete points and x = 0 will have its temperature stored in T[N];

float T[MAX], Tn[MAX];  
 
...  
 
 for(t=1;t<steps;t++){  
update(T,Tn,coeff1, coeff2,MAX);  
swap(Tn,T,MAX);  
    }  
...  
 
void update ( float *Told, float *Tnew, float c1, float c2,int size){  
 int i,M;  
    M=(size-1)/2;  
    for(i=-M+1;i<=M-1;i++){  
*(Tnew+i+M)=c1*(*(Told+i+M))+c2*( *(Told+i+M-1) + *(Told+i+M+1));  
    }  
/* now get the ends */  
*(Tnew+0)=c1*( *(Told+0)) + c2*( *(Told+1));  
*(Tnew+size-1)=c1*( *(Told+size-1)) + c2*( *(Told+size-2));  
}  
 
 
void swap ( float *Told, float *Tnew, int size){  
    int i,M;  
    M=(size-1)/2;  
    for(i=-M;i<=M;i++){  
*(Tnew+i+M)=*(Told+i+M);  
    }  
 
}  

Graphics programming with GNU libplot

In order to use the library one first compiles and installs it, and then includes its header file in the program preamble

#include<stdlib.h>  
#include<math.h>  
#include<plot.h>

and link with the library, as well as Xt, X11, ICE, Xext and several other Xwindows libraries

gcc stepwave.c  -lplot -lXaw -lXmu -lXt -lSM -lICE -lXext -lX11 -lm

Plotter handle (file descriptor)

One of the variables that needs to be declared in order to use the library is a file descriptor for the plotter device, be it GIF, X, PS or one of the other choices. We will create animated pseudo-GIF89a movies, output to the Xplotter or create PostScript graphics with the three devices listed above, or else use the PNM device to create pnm files suitable for encoding into an MPEG stream. At any rate the preamble of our program now looks like

#include<stdio.h>  
#include<strings.h>  
#include<stdlib.h>  
#include<math.h>  
#include<plot.h>  
 
int pl_handle; /* plotter file descriptor */  
 
double firstX, firstY, secondX, secondY,  
int n,m,l,time;

Choosing the plotter device with pl_newpl()

Libplot is extremely sophisticated in how it sets up a plotting window using affine or vectorized graphics in a user-defined coordinate system. We set up the plotting device by specifying its parameters such as size, type, backround color, transparency, and so forth, and then open the plotter by initializing its file descriptor with pl_newpl(). The first argument passed to it is the plotter device, “X”, “gif”, “ps”, “pnm” are some of your options.
For example we take the following steps to open an animated GIF plotter with a 5 ms frame delay and transparent white background. This will make any white objects drawn in the plotter transparent, so that when viewed in a web browser, which knows how to load GIF89 images, the background does not appear.

 fptr = fopen ("well.gif", "w+");      /* GIF written to this file */  
  pl_parampl ("BITMAPSIZE", "300x150");  
  pl_parampl ("BG_COLOR", "white");  
  pl_parampl ("TRANSPARENT_COLOR", "white");  
  pl_parampl ("GIF_ITERATIONS", "100");  
  pl_parampl ("GIF_DELAY", "5");  
  pl_handle = pl_newpl ("gif", stdin, fptr, stderr); /* actual initialize */  
  pl_selectpl (pl_handle);            /* select this plotter for output */  
  if (pl_openpl () < 0)         /* open Plotter */  
    {  
      fprintf (stderr, "Couldn’t open Plotter\n");  
      return 1;  
    }

All of the plotters parameters must be set with pl_parampl() commands before the file descriptor is declared (the plotter device is opened).
The code block above opens a gif plotter with our desired parameters, with output into file fptr and error output to stderr, which is the console. We could declare several plotters in the same program and switch between them with the pl_selectpl() command. To open the Xdisplay instead and double buffer, or build the next image while displaying its predecessor, with

  pl_parampl ("BITMAPSIZE", "300x150");  
  pl_parampl ("BG_COLOR", "white");  
  pl_parampl ("USE_DOUBLE_BUFFERING", "yes");  
  pl_handle = pl_newpl ("X", stdin, fptr, stderr);  
  pl_selectpl (pl_handle);  
  if (pl_openpl () < 0)         /* open Plotter */  
    {  
      fprintf (stderr, "Couldn’t open Plotter\n");  
      return 1;  
    }  

This is the initialization that we would use to perform a real-time animation in Xwindows. For truly numerical-intensive simulations, real-time animation could be painfully slow, and in those cases we use the GIF89a devices, or output snapshots of the action as PNM graphics to be encoded into an MPEG stream.

The plotter window geometry

Once the plotter has been initialized we set up an affine coordinate system within it with the command pl_fspace(). This takes four floating point arguments, the (x,y) coordinates of the lower left corner of the plotting window, and the (x,y) coordinates of the upper right corner. For example, suppose that we have created and stored a large set of data points ı

(xn,yn) = (point[n][0],point[n][1])
and we wish to create an affine space large enough to plot all of them with some room to spare, we would set up the plotter as follows
printf("\tDetermining the scale of the snapshots;\t");  
 
largest=point[0][0];  
 for(n=0;n<N;n++){  
if(point[n][0]>largest)largest=point[n][0];  
if(point[n][1]>largest)largest=point[n][1];}  
 
 pl_fspace (-1.1*largest,-1.1*largest,1.1*largest,1.1*largest);

There is no need to take into account aspect ratio if the X and Y-limits of the window are of the same size.

Graphics primitives and attributes

We now specify various graphics primitive attributes such as linewidths, text sizes, object colors and fill types. The following example illustrates a variety of these attributes and the commands used to set them

 pl_filltype (1);                   /* fill in the object */  
 pl_pencolorname ("red");  
 pl_fillcolorname ("red");  
 pl_fcircle(0.0,0.0,0.04*largest);   /* draw the nucleus      */  
                                     /* centered at (0,0)     */  
                                     /* radius is 4 percent   */  
                                     /* of coordinate window  */  
 
 pl_fmove(-0.9*largest,-0.9*largest); /* move to this location in window */  
 pl_ffontsize(0.1*largest);           /* textheight 1/20 window width    */  
 pl_label("3,2,0");                   /* and draw this label string      */  
 pl_filltype (0);                     /* subsequent objects not filled   */  

The actual routines to create graphics primitives all begin with pl_, and come in two flavors, we want to use those that begin with pl_f, which will place objects at locations whose coordinates are floating point numbers, rather than integers. We can create lines, points, circles, ellipses, boxes, bezier curves, labels and a variety of other objects. We can specify colors with commands such as

pl_pencolorname("gray60");  
pl_pencolor(55000, 65000,34000);  
pl_fillcolor(0,0,65535);  
pl_filltype(1); /* fill=yes */  
pl_linemod("dotted");  
pl_linemod("solid");  
pl_fillcolorname("yellow");  
pl_ffontsize(0.35);  
pl_fontname("HersheySerif");

and so forth. For the full list of commands and syntax, read the info page

less /usr/info/plotutils.info  
or  
less /usr/local/info/plotutils.info

Closing the plotter

The final act in a libplot program is to close the plotter device and de-allocate its file descriptor with

  pl_closepl ();  
  pl_selectpl (0);  
  pl_deletepl (pl_handle);

You will find it extremely easy to program high quality graphics and animations in X with this library. I recommend the 2.2 version plotutils-2.2.tar.gz over other revisions.

Make and the Makefile

Make is a utility for automating the building and re-building of large programming projects that contain many sources. If you have placed all of your source-files ( *.c files) in a directory, and type “make” in that directory, make will look for a special file, called the Makefile, that consists of substitutions or identifications, targets, dependencies and directives, and execute the instructions in the file to build the project.
The advantages to using make are many; you can define long compiler invocations in the Makefile rather than type them in at compile time, you can automate compilation of a project that contains thousands of C source files, and if you change one or two of your sources, make will re-compile only those sources that depend on the files that you changed, rather than re-compile the entire project.

Here is the Makefile for the four versions of our heatflow project;

# Simple C Makefile for 290 projects  
CC = gcc  
CPP = /lib/cpp  
 
FLAGS = -Wall -fomit-frame-pointer -O2  
LDFLAGS = -I./ -L./ -lm  
XFLAGS = -L/usr/X11R6/lib -lplot -lXaw -lXmu -lXt -lSM  -lICE -lXext -lX11 -lm  
 
# targets  
all:    heatflow1 heatflow2 heatflow3 heatflow4 \  
 
heatflow1: heatflow1.c  
        $(CC) $(FLAGS) heatflow1.c $(LDFLAGS) -o heatflow1  
heatflow2: heatflow2.c  
        $(CC) $(FLAGS) heatflow2.c $(LDFLAGS) -o heatflow2  
heatflow3: heatflow3.c  
        $(CC) $(FLAGS) heatflow3.c $(XFLAGS) -o heatflow3  
heatflow4: heatflow4.c  
        $(CC) $(FLAGS) heatflow4.c $(XFLAGS) -o heatflow4  
 
clean:  
        -rm  heatflow1 heatflow2 heatflow3 heatflow4

Make is very sensitive about syntax; the dependency-directive pairs have the directive lines beginning with a tab, not eight spaces. If you get a “missing seperator” error from make, check your file for tabs. I find that the vi editor is the best one for writing Makefiles.

To build all of the targeted executables, type “make” or “make all”, to build just one target, say heatflow2, type “make heatflow2”, and to remove stale compiled targets for a fresh build, type “make clean”.

Lines that begin with # are comments, and are ignored. Values of defined objects or substitutions are referred to by prepending the substitution with a dollar sign, just like you would do in a shell script.