VI.1.1 Introduction

Mathematica is an excellent environment in which to create simulation programs. The combination of a superb high-level programming language, an interpreted environment that encourages step-by-step development and testing, and rich visualization tools provides unparalleled ease of development. The only drawback is that Mathematica is demanding on system resources, both speed and memory. Fortunately, Mathematica has a facility called MathLink that makes it easy to integrate external functions written in a compiled language (like C, for example). With MathLink it is possible to write external programs that become seamless extensions to the Mathematica kernel. In this way, if you have an algorithm that needs to be written in C for efficiency reasons, you still have access to the rich Mathematica environment for all the auxiliary tasks -- prototyping, preparing and managing input, analyzing and visualizing output, etc.

MathLink works in two directions: you can use existing C functions in a Mathematica program, and you can use existing Mathematica functions in a C program. Calling C functions from within a Mathematica program cannot only improve the overall performance of a simulation, but also may provide the user with the means to access external code not directly accessible from Mathematica, like numerical libraries. From this point of view the C code acts as a back-end for Mathematica. The other direction, in which Mathematica functions are called from within a C program, can be very useful when a user wants to utilize the visualization capabilities of Mathematica. The C program then generates the data to be fed into Mathematica to display the program's results. From this point of view Mathematica acts as a presentation tool of the C program.

This notebook will only discuss one of the two above-mentioned directions. For the sake of simplicity of explanation we will describe how to use MathLink to call external functions from Mathematica. After a short basic introduction to MathLink, a very simple step-by-step example will illustrate how to connect a Mathematica program with a C program. The notebook will be concluded with a more meaningful application: the physical problem of heat transfer in a rod.

We assume a basic familiarity with C, and since our focus is on using MathLink we will not deal extensively with the actual implementations of the algorithms, except to highlight MathLink programming issues and contrast programming styles in Mathematica and C.

The source code is ANSI C, and it will compile without modification on any platform Mathematica supports. Instructions for building and running the programs appear below.

The code presented here compiles under version 2.2.2 of MathLink (the current shipping version on all platforms as of this writing) and later.

It is not our intention here to provide an in-depth tutorial on MathLink. The main documentation for MathLink is the MathLink Reference Guide, which ships with most versions of Mathematica and can also be obtained on MathSource or ordered from Wolfram Research directly. In addition, there is an extensive MathLink tutorial on MathSource, as well as various sample programs. You should consult these materials for more information about MathLink programming.

MathLink is a library of functions that implement a protocol for sending and receiving Mathematica expressions. The easiest and most common application is to allow external functions written in other languages to be called from within the Mathematica environment. This is how we use MathLink in writing the simulations. As we shall see, it is a relatively simple matter to incorporate such routines into Mathematica.

We refer to external functions that are called from Mathematica as "installable" functions, since they use the Install function from Mathematica to be incorporated into the Mathematica environment. For each external function you want to call from Mathematica, you write a "template" that specifies information about the function, including its name, the arguments that need to be passed and their types, and the type of argument it returns. This template file is then passed through a tool called mprep, which writes C code that manages most of the MathLink-related aspects of the program. The magic here is in this processing by mprep. It writes a lot of C code that "wraps around" your function and handles a lot of complexity you would rather not have to deal with. It also lets you write portable programs, because all the machine-specific code is written by mprep according to what platform it is running on.

VI.1.2 A simple example

Suppose we want to call from Mathematica a C function that returns the sum of two integers. The C program below would do the job if Mathematica would not be involved.

    int main (int  argc,  char *argv[]) {
        int p;
        p = add (10, 20);
        return 0;
    }
    
    int add (int  m, int  n) {
        return (m + n);
    }    

When we want to call the function add from Mathematica the above C program has to be modified slightly. The main function must be replaced, and each function that Mathematica is to call must be annotated with a MathLink template.
Because Mathematica will be calling add directly, there is no more need for the original version of main. Instead, it will be replaced with some lines required for the proper use of MathLink. The new and simple version of main is

    int main (int  argc,  char *argv[]) {
        return (MLMain (argc, argv));
    }

The only function from the original program that is to be called from Mathematica is the add function. The annotation of this function does not change the function itself. Rather, it proceeds a template to the function needed to tell Mathematica everything that it needs to know in order to call add. The template contains the function's name, its return type, and the name and type of each of its parameters. Though some of the information in the template is redundant, all of it is required. The annotated version of add looks like

    :Begin:
    :Function:        add
    :Pattern:        add[m_Integer, n_Integer]
    :Arguments:        {m, n}
    :ArgumentTypes:    {Integer, Integer}
    :ReturnType:    Integer
    :End:
    
    int add (int  m, int  n) {
        return (m + n);
    }

Note that the (Mathematica) type name Integer stands for the (C) type int. Likewise, the type name Real stands for double etc.

To complete the annotated program we should include the mathlink.h header file for a number of (type) definitions and macros. The complete program with the annotations for MathLink is given by

#include <mathlink.h>

int main (int  argc,  char *argv[]) {
    return (MLMain (argc, argv));
}
    
:Begin:
:Function:        add
:Pattern:        add[m_Integer, n_Integer]
:Arguments:        {m, n}
:ArgumentTypes:    {Integer, Integer}
:ReturnType:    Integer
:End:
    
int add (int  m, int  n) {
    return (m + n);
}

Annotated version of C program

The annotated program is no longer a regular C program due to the presence of the function template. To distinguish the extended C code from a normal C program, it is custom to put the annotated code in a file with a postfix .tm (think of "tm" as standing for template). The program above could for example be placed in a file called simple.tm.
The next step is to translate the annotated program into valid C code. For this purpose the Mathematica environment provides the mprep translator. The result of applying the translator, a valid C program, can then be compiled to create an executable. For convenience, Mathematica also presents the user with the mcc script, which includes both the translation and the compilation. To compile the simple.tm file and create an executable named simple we just give the command

    mcc simple.tm -o simple -lsocket -lnsl

at the prompt in the directory where the simple.tm file resides.
Once we have created the executable the last step to take is to install it into Mathematica. Within Mathematica the executable is installed by the command

Mathematica is now able to call the add function just as if it was part of Mathematica itself. The following expressions are valid

VI.1.3 Heat transfer in a rod

A more serious example of the use of MathLink is discussed in this subsection. Consider a rod with a constant cross section made of some heat conducting material, for example silver or copper. The rod is completely wrapped up in some thermal insulator except for the end points, and is initially kept at constant temperature. In this situation the only points where heat can enter or leave the rod are its end points. At some moment heat sources of different temperature are connected to both end points. If we assume that a heat source of 100 °C is attached to the left end of the rod and a heat source of 50 °C to the right end, the rod will start warming up, more quickly at its left end then at its right end. However, the temperatures at the end points of the rod remain constant. After some time we will have a situation in which the temperature in the rod linearly varies from 100 °C on the left to 50 °C on the right end. The temperature of the middle point is then 75 °C. What is really interesting to know is how the temperature in the rod changes as a function of time, for example in the middle point.

In the context of heat transfer problems the name of the French mathematician and physicist Jean Baptiste Joseph Fourier (1768 - 1830) should be mentioned who devoted a substantial part of his life to the study of the physics and mathematics of heat transfer.

Before a model can be developed for the heat transfer in the rod a few assumptions have to be made. We have already assumed that the rod is encased in thermal insulation and that its initial temperature is constant. Also, we have assumed that the cross section was the same throughout the rod. In addition one extra assumption is made: the rod is made of homogeneous material, i.e. is made of one material with only one set of heat conducting properties.

Within the model to be derived a lot of properties and parameters are involved. The length and cross section of the rod, the initial temperature, the thermal conductivity (how much heat is transferred per degree Celsius in one length unit of the rod during one second), the density, the specific heat (how much energy is required to heat one mass unit of the rod one degree Celsius) and the thermal diffusity (the area through which heat spreads during one second) are all relevant for a realistic modelling of the rod. The temperatures at the sources also play their role within the model. In the frame below a glossary of all parameters used in the model are listed.

Rod characteristics

L    length (cm)
A    cross-sectional area ()
T    initial temperature
κ    thermal conductivity (cal / sec / cm / °C)
ρ    density (gm / )
σ    specific heat (cal / gm / °C)
    thermal diffusity (κ / (ρ σ) / sec)

Heat Sources

    temperature of left source (°C)
    temperature of right source (°C)

Discretization

n    number of segments into which rod is divided
h    L / n, length of each segment (cm)
Δt    discrete simulated time interval (sec)
    temperature at segment boundary i = 0 ... n at time t (°C)
    temperature at segment boundary i = 0 ... n at time t + Δt (°C)

Model parameters

The model that will be developed shortly is based on three basic facts from theory. Two of these say something on what will happen once the rod reaches steady state, whereas the third says something on the situation before steady state. The facts are:

1.  If we leave the rod in contact with the heat sources, it will eventually reach a steady state. At steady state, the temperature along the rod will vary linearly between the temperatures of the two heat sources. This means that the temperature halfway between the two end points will °C.  (VI.1)

2.  Once the rod reaches steady state, heat energy will flow from the hotter heat source, through the rod, to the cooler heat source. The amount of heat that will flow from left to right during an interval of Δt seconds will be
                     calories.    (VI.2)
    This is known as Fourier's law of conduction. The quantity in the expression is the temperature difference per unit length and is called the temperature gradient.
   
3. When heat enters the rod, the rod's temperature will rise. If the energy entering is ΔH calories, the rod's temperature will rise by
                         °C.    (VI.3)
    Note that
  Note that ΔH is a signed quantity, so if heat leaves the rod, its temperature will fall.

The model

The approach to tackle the heat transfer problem is by means of discretization. We start with dividing the rod into n separate segments, each of length . The figure below shows the case where and each segment has length .

The six segments form seven boundaries: five internal boundaries and the ends of the rod. The temperatures at the boundaries are labelled through , where is the same as (the temperature at the rod's left end) and is the same as (the temperature at the rod's right end). We will assume that each of the segments is at steady state. The assumption that each of the segments is at steady state may seem a little unrealistic at first sight. However, by dividing the rod into a large number of segments the errors introduced by the assumption can be minimized. In fact, it is even less realistic to think that the entire rod is at steady state!
Under this assumption the temperature in each of the segments varies linearly between the temperature at its left end to that of its right end, e.g. in the first segment the temperature varies linearly between and . In the general case with segments, at the temperatures at the segment boundaries are given by and , for . In each time step heat will flow among the segments and the temperatures at the interior boundaries will change. By assuming that the segments are in steady state, we are in a position to use the equations discussed in the previous subsection to model the heat flow between the segments and to determine how the temperatures at the boundaries change.
If the temperatures at the segment boundaries at some time are given by through , what will be the temperatures , ..., at these boundaries at for some small Δt, under the assumption that the individual segments are at steady state? The fact that the rod's end points are connected to the heat sources implies that the temperatures and will not change in any time interval. The points of interest are therefore the interior segment boundaries. In the following figure a detailed view is given of two contiguous segments of the rod.

The temperatures at the three boundaries are , and respectively. Each of the segments have been divided in half (the dashed lines), and each of the resulting half-segments has been labeled (a through d). As the lengths of the segments is simply , the lengths of the half-segments is simply given by .
By determining how much heat flows into the half-segments b and c during the interval , we can derive how the average temperature changes within these two half-segments and consequently at their shared boundary. What we know is that heat will flow across a into (or out of) b, and across d into (or out of) c. From equation VI.1 (see previous subsection) it follows that the temperature at the boundary between a and b is given by . From equation VI.2 it easily follows that the amount of heat flowing across a into b during the interval Δt will be

                 .

Similarly, the heat flowing across d into c will be

                 .

Combining these two results, the total heat change in the half-segments b and c will be

                 .

Because the total length of the two half-segments b and c is equal to h, equation VI.3 tells us that the change in the average temperature of the two segments will be

                 .

From this equation an expression for can be derived immediately: .

Recapitulating, we find that

	=
	=
	=

where is called the thermal diffusity of the rod material.

Numerical Algorithm

The model that was described in the previous subsection is a so-called finite element model. The essence of the model was the division of the rod into separate elements (segments) followed by applying physical analysis to the elements. In addition, a simplifying assumption was made on the steady state of the segments. Applying this assumption to the entire rod would have been inaccurate, but for the segments it is valid as long as the lengths of the segments are small. For small segments the errors introduced will be acceptable and by increasing the number of segments the amount of error decreases. In the limiting case with an infinite number of segments the exact solution can be reached.
It may be obvious that in a program the number of segments is finite. The expressions at the end of the previous subsection enable us to determine the temperatures at the segment boundaries at time once the temperatures at time are known. For a given initial temperature of the rod it is thus possible to calculate the temperatures at the segment boundaries at any time (with resolution ) by repeatedly applying the above equations.
Examining the equations we see that the following quantities need to be known in our model: the length of the rod, the (constant) temperatures at its end points, the initial temperature of the rod itself, the thermal diffusity of the material, the number of segments into which the rod should be divided, and the length of the time step Δt. By choosing a convenient value for Δt, we can even simplify the equations. If we choose Δt such that

the last equation is reduced to

                 ,     .

The C program, simulating the heat transfer in a rod, is shown below. The main function does nothing more than asking for the model's parameters followed by calling the function simulate. In the version below, this function returns the approximate temperature in the middle of the rod after the specified simulation interval has elapsed. Notice that in this implementation the maximum number of segments is limited to 1000. To increase this number the user must modify the MAXSIZE value before compilation.

#include <stdio.h>

#define MAXSIZE 1001    /* The maximum number of segments that can be    */
                        /* simulated.  This is all that must be modified */
                           /* to change the maximum size of the simulation. */


/* This function models heat transfer in a rod of uniform temperature
* that is insulated along its lengths and has a heat source applied
* to its ends.  The rod's thermal "diffusivity" in cm^2/sec and
* "length" in cm are given, as are the temperature of the "rod", the
* "left" heat source, and the "right" heat source in degrees Celsius.
*
* The function returns the approximate temperature of the center of
* the rod in degrees Celsius "duration" seconds after the heat sources
* are applied.
*
* The implementation is based on a finite-element model with "segs"
* elements.
*/

double simulate (double diffusivity, double length,
                 double left, double right, double rod,
                 double duration, int segs)
{
  double delta;         /* Length in seconds of each simulated time step */

  double cur[MAXSIZE],  /* Temperatures at the ends of each of the "segs" */
                        /* segments, where cur[0] is the temperature of */
                        /* the left heat source and cur[segs] is the */
                        /* temperature of the right heat source. */

         nxt[MAXSIZE];  /* Used to temporarily hold the new temperatures */
                        /* of "cur" during the simulation. */

  int i;                /* Steps through "cur" and "nxt". */

  /* Compute the length of the simulated time step. */
  delta = (length*length/(segs*segs)) / (2*diffusivity);

  /* Initialize "cur" to reflect the temperatures at time zero. */
  cur[0] = left;
  for (i = 1; i < segs; i = i+1) {
    cur[i] = rod;
  }
  cur[segs] = right;

  /* Simulate the passage of time.  During each iteration, compute the
     new segment boundary temperatures in terms of the old segment
     boundary temperatures. */
  while (duration > 0) {
    for (i = 1; i < segs; i = i+1) {
      nxt[i] = (cur[i-1] + cur[i+1]) / 2;
    }

    for (i = 1; i < segs; i = i+1) {
      cur[i] = nxt[i];
    }

    duration = duration - delta;
  }

  /* Return the temperature of the center segment boundary. */
  return(cur[segs/2]);
}


/* Main program */

void main (void)
{
  double C, length, left, right,  /* Parameters specified by the user. */
         rod, time, temp;
  int segs;                       /* Size of finite-element simulation. */

  /* Read constants from the keyboard using "scanf". */

  printf("Enter thermal diffusivity (cm^2/sec): ");
  scanf("%lf", &C);

  printf("Enter length of rod (cm): ");
  scanf("%lf", &length);

  printf("Enter temperature at left end (degrees C): ");
  scanf("%lf", &left);

  printf("Enter temperature of rod (degrees C): ");
  scanf("%lf", &rod);

  printf("Enter temperature at right end (degrees C): ");
  scanf("%lf", &right);

  printf("Enter length of simulation (sec): ");
  scanf("%lf", &time);

  printf("Enter number of segments: ");
  scanf("%d", &segs);

  /* Calculate and display the result. */

  temp = simulate(C, length, left, right, rod, time, segs);
  printf("Temperature at center after %g seconds is %g degrees C\n",
         time, temp);

}

C implementation of the head transfer in a rod

The function simulate reflects the model as discussed above directly.  After setting and initializing the begin temperatures for the rod and its end points, the temperatures for the segment boundaries are computed over the simulation interval. Finally, the end temperature for the middle of the rod is returned.

REQUIRED
VI.1.a Run the above C program. Take a silver (thermal diffusity: 1.752 ) rod of 10 cm., which is initially at 0 °C. Let the left and right heat sources have a temperature of 100 °C and 50 °C respectively. How long does it take for the rod's middle point to reach its end temperature?

REQUIRED
VI.1.b Modify the C program in such a way that the function simulate can be called directly from within Mathematica. (See the simple example in a previous subsection on how to use Mathematica and C together.)

REQUIRED
VI.1.c Modify the function simulate in such a way that it returns a set of values with the temperatures at the middle of the rod for each time step, i.e. the returned set should contain the temperatures at for the middle of the rod. The return type should be the array of floats (reals). Plot in Mathematica the temperature at the middle point versus the time.

REQUIRED
VI.1.d  Adapt the program from VI.1.c so that the function simulate returns the temperatures at    for all segment boundaries. Plot the result for a single time step in Mathematica using a color map, where the rod is drawn in a 2-D figure and the temperatures at the segment boundaries are shown with colours ranging from blue (low temperature) to red (high temperature). Do this for a number of time steps and generate an animation of the heat transfer in the rod.