Documentation - HMD Mathematics

Login as the root user, or su root. Copy the HMD package file to the /usr/local/src directory.

Change directory to /usr/local/src and unpack the HMD package there.

Replace "X.X.X" shown above with the version number, e.g., 4.1.2. A new directory will be created under /usr/local/src and the files will be copied there.

Change directory to the new directory just created.

Build the project by entering the following commands.

After the "make install" step, but only after the "make install" step, type

Leave the project in /usr/local/src, because you might want to uninstall it later. To uninstall, change directory back to the hmd-X.X.X directory and type "make uninstall".

An interactive session can be started with the "-i" command line flag.

Or the HMD program can be run in batch mode in which a script file is executed, after which the program exits.

Script commands execute statements or function calls. Statements are either mathematical assignments in the form of an equation, or are simple commands. Her is an example of making a simple assignment.

A variable named x is created and assigned the value of 1. A variable created this way, that is, without declaration, will automatically be given the type of double precision.

Now a mathematical calculation can be performed,

Arrays

A variable can be created as an array,

And individual elements of the array can be changed,

or

We can access a range of array elements,

We can put a range of indices into a variable and use that variable as an index into an array,

Arrays can be built by combining other arrays,

Notice that the comma character separates array elements in the same row. The semicolon within an array separates multiple rows.

Enumeration Arrays (or range indices)

Arrays can be defined by an enumeration within a range of values. The enumeration is defined by giving a starting number, a colon separator, and an ending number.

A two-point enumeration as defined above is range of values starting with the first number, where each subsequent array element is produced by adding 1. The last element of the enumeration will be the greatest number less than or equal to the ending specifier.

An enumeration array with a fractional subdivision can be created using three specifiers separated by colon characters. The syntax is start : increment : end.

The increment is allowed to be a negative number, which allows counting backward.

Incrmentally Building Arrays

It is often convenient to build an array, as a literal object, in a loop by starting with an "empty array." Then, with each iteration of the loop, append another element of the array. An empty array is first created using the nil token in square brackets.

Data Types

Variables may be declared, as in the C langauge. Or they may be defined instantly by assigning a value to a new symbol. The numerical data types can be declared in the HMD workspace:

• char

  An 8-bit signed character, like in C.

• uchar

  An 8-bit unsigned character, like in C.

• short

  A 16-bit signed integer, like in C.

• ushort

  A 16-bit unsigned integer, like in C.

• int

  A 32-bit signed integer, like in C.

• uint

  A 32-bit unsigned integer, like in C.

• long

  A 32-bit signed integer, like in C.

• ulong

  A 32-bit unsigned integer, like in C.

• float

  A 32-bit floating-point, like in C.

• double

  A 64-bit floating-point, like in C.

• complex

  A complex number composed of a 32-bit real part and a 32-bit imaginary part.

• doublecomplex

  A complex number composed of a 64-bit real part and a 64-bit imaginary part.

The syntax for declaring a variable resembles a declaration in the C language. For example, to declare an integer the following code would be used.

Note that if a new variable is created by an assignment to a new symbol its data type will be double, for single numbers as well as arrays.

Assignment to a complex type, complex or doublecomplex is performed by enclosing two numbers in angle brackets separated by a percent-i symbol.

A string is created either with a declaration or by just assigning a variable.

Redirection - IF / ELSE

Redirection of script execution is done with the if / else block.

Several rules must be followed for an if / else block.

• The condition to be tested is enclosed in parentheses.
• The body of the redirected code is enclosed in curly brackets.
• The opening curly bracket must be on the same line as the if /else keyword.
• The closing curly bracket must be on a line by itself.
• There is no else if combination, like in C. The same thing must be accomplished with an if inside an else block.
• Redirection is not possible from the interactive command prompt. In other words, multiple-line code blocks cannot be entered. These are executed from within script files.

The else block, when needed, must be given on the next line,

New variables created within an if / else block will become global variables. However, a local variable may be declared inside the code block, in which case it will not exist globally. It will cause an error to try to access such a variable outside the code block.

Redirection - WHILE loop

Looping a section of script commands is accomplished with the while code block. This resembles the C language while loop.

Like the if / else blocks, a while block may contain local variable declarations.

Observe that the while loop must conform to similar syntax rules as the if / else block.

Redirection - FOR loop

In addition to the while loop, a for loop redirection is provided, which is sometimes more convenient than a while loop. The counter variable in a for loop does not need to be declared or initialized outside the loop.

This for loop construct is different than other programming language syntax in that the test condition is not a comparison operator. It is, rather, a list containing the loop variable name and the initial and final values. Also, the values in parentheses are separated by the DOT-COMMA operator. In the above example the first entry, j, is the loop counter. The second entry, 1, is the initial value for j, and the third entry, 3, is the final value of j. The loop counter, j, becomes a local variable within the loop block.

Unlike other for loop constructs there is no choice about the comparison attribute on the final value of the loop counter. The logic is always "less than or equal," that is, the loop continues while j is less than or equal 3. Also, be careful that there is always whitespace between the DOT-COMMA operators and the symbols they separate.

An optional fourth, increment, value may be specified as the last element in the for list.

The for loop conforms to the following rules of syntax:

1. The first list element is a variable name for a local variable that will serve as the loop counter inside the loop block.
2. The second list element must be the initial value for the loop counter. This may be either a literal number or an already-existing variable containing a real-valued number. If the initial value is a variable it will only be accessed before the first iteration through the loop.
3. The third element in the list must be the final value of the loop counter. The loop will continue to iterate as long as the counter variable is less than or equal to the final value, when counting up. When counting down the loop will continue as long as the loop counter is greater than or equal to the final value. The final value may be either a literal number or a variable. If the final value is a variable then its value can be changed inside the loop.
4. There can be an optional fourth list element which is the increment value. The increment value may be either a literal number or a variable. If the increment is given by a variable then its value can be changed inside the loop. When an increment value is given, then its sign determines whether the loop counts up or down. If no loop increment value is given, then the loop increments by whole numbers by default, and the direction of the iteration, up or down, is determined by the magnitudes of the initial and final values.

Many scripting languages have a capability that could be described as a "foreach" loop, in which a loop counter sequentially reads the values of an array or a list. This script interpreter can do that by using the in operator in the for conditional clause.

The for loop can be used with compound conditional statements that create a single logical AND as a condition. That way several loop counters can be iterated simultaneously.

The AND operator && forces the loop to exit when the first for loop condition becomes false. This will prevent the loop variable that iterates over an array from exceeding the index bounds of the array. If you put an OR condition, that is, || as a compound loop connective the loop will execute, but the behavior is undefined.

Command History

At the interpreter command line (parser command entry) the command history will be saved for the next session if the Gnu Readline library was included in the build. If the Gnu Readline package was already installed when you build HMD, then the HMD build will include the Readline library by default. If not, then install Readline and uninstall, reinstall HMD.

The command history is saved to a local file in the current working directory. The file name will be ".hmd_history".

Help

Command line help is available by typing

A list of subtopics will be displayed. Any subtopic can then be displayed by typing help and the subtopic in single quotes.

Functions can be defined in a script file and executed later, either at the command prompt or in another script file. Here is a simple function definition.

When working from the interactive prompt, the above definition cannot be typed line by line into the interpreter. It must be contained in a script file, then read into the workspace with the source command:

A function definition must begin with the function keyword, followed by the function name, followed by an opening curly bracket, followed by zero or more input parameters, followed by a closing curly bracket, followed by whitespace and the opening curly bracket for the function body. These all must be on the same line.

Functions are defined for the purpose of providing relocatable code blocks within the script. The function parameter names will be created in an independent namespace than global parameters. Inside the body of a function local variables may be declared, and these variables will be created in a separate namespace than the global variables. These local variables will not exist outside the function body.

In this example the following script is written into a file "my_func.hmd" and then run as a batch file:

Now running the script from the unix prompt:

Functions are allowed to be recursive. The following script is a factorial function.

And running it as a batch file,

To force the scope of a variable to be local within the function, the variable must be locally created with a declaration. Any variable created by assignment will be global.

The primary file type is an HMD script file, which can either be read and executed as a batch file from the Unix command line or can be executed within the workspace interpreter using the source command.

Data files must contain rectangular arrays and matrices of real or complex numbers. For example, the following matrix of complex numbers can be read into a workspace matrix with the read function,

The following script line reads the matrix, contained in a file called "read_keep2.dat" into a matrix variable named X.

The numbers in each line (row) of the data file must be separated by a space, tab, or comma. If the numbers are complex, then there cannot be any whitespace between the + symbol or the i symbol.

The read function can take an optional second argument: either a data type name, or the string 'block'. The data type names can be one of any of the following: 'char', 'uchar', 'short', 'ushort', 'int', 'uint', 'long', 'ulong', 'float', 'double', 'complex', 'doublecomplex'. If not data type specifier is given then the default data type will be either double or doublecomplex, depending on the data in the file.

Example,

If the second argument to the read function is the string 'block', then the data will be read into a special matrix variable (called a block) that is internally simpler than a standard matrix or array. A block may be only of type double or doublecomplex, and may be used only for a limited number of operations.

For saving a matrix to a file there is a corresponding function called save.

When saving a matrix the numbers in the data file will be written in general numerical format, meaning that the least number of significant digits will be shown. It is sometimes desirable to print numbers with a fixed width, which makes columns align. In that case the default file save format can be changed with the save_matrix_fmt_g global variable. The default value is 1, meaning that general format will be used. To use fixed width (of 1.4e) set this variable to 0.

HMD has a plot function which is a front end to the Gnuplot program. Gnuplot is provided as a standard software package on most Linux distributions. It has become the defacto universal plotting program in the Linux world. For the plot function to work the Gnuplot program must be accessible in the PATH. You can test this by opening a command terminal window and typing gnuplot.

To plot a single data array as an X-Y plot just put the name of the data array as a single argument to plot. For example,

Or, you may want the x data values displayed on the x axis, not just the data indices. Then, the x data set is provided as the first argument followed by the y data argument.

You can plot more than one data set on the same plot. You must give the x data as the first argument followed by multiple y data arguments.

You can change the style of any one of the y data curves by supplying an options argument after the y data argument. If there are multiple y data arguments then an option string may be given after each y data array. The syntax is ydata, options, ydata, options, ...

The options for a curve (a y data array) are provided as a string, which is delimited by single quotes. The options format follows exactly that which would be provide if running the Gnuplot program by itself. Further documentation of the options can be found in the Gnuplot documentation and will not be duplicated here.

Gnuplot is a powerful program which allows you to set up a custom configuration for your plot. These are option commands that would be given to Gnuplot before the plot command is executed, either on the command line or from a script file. The HMD plot function provides a way to give set up options for the plot as a single string which must be the first argument. So the complete syntax for the plot command is

plot([setup], x_data, y_data_1 [, options_1], y_data_2 [, options_2] ...)

In the following example a set up option is given to turn on grid lines.

Tp provide more than one setup option these must be given within the single setup string separated by semicolons.

Introduction

The term "finite difference" is sometimes used to describe approximation techniques for solving differential equations that involve stepping a difference equation. But the subject of finite differences has traditionally been understood to mean the fundamental principles of numerical mathematics, if not the foundations of The Calculus. As such, there are certain mathematical procedures that are expected to be included in the subject, such as Difference Tables, the Difference Calculus, the Summation Calculus, Stirling Numbers, and Difference Equations.

The HMD finite difference library (generally found in finidiff.c) is intended to provide some of these mathematical techniques as script functions.

Difference Tables

The difference table can be used for studying the properties of a data set, such as evaluating roundoff error, as well as constructing interpolation functions. The HMD function difftab_fndf takes an array of data, calculates the difference table, and returns a handle to an internal object that stores the table.

For example, the following code creates and prints the table,

The print_difftab_fndf provides a convenience for displaying the difference table in the usual format with differences lying between the previous values in the table

or, if an addition argument of "1" is given to the function, the table will be oriented sideways,

The quality of your data samples, that is, the smoothness, or even the continuity, may be studied with its difference table. The capability of the difference table to do that is related to The Fundamental Theorem of the Difference Calculus:

It is related to the fact that the (N+1)th derivative of an Nth-degree polynomial is zero, which has far-reaching implications for the concept of a continuous-valued function.

This means that for a continous-valued function each subsequent set of differences in the difference table should have absolute values that become smaller than the previous set. When subsequent differences become larger instead of smaller it is similar to high frequency components in a Fourier series, which can reveal a discontinuity in the function. For numerical mathematics the discontinuities are usually indications of errors, such as roundoff, in the data.

In the case of roundoff error it can be shown (Hamming, section 10.3, "Error Propagation in a Difference Table") that roundoff error in a data set will result in alternating signs in the fourth or fifth difference in its difference table. For example, consider a difference table of random numbers (which is 100 percent error),

the difference valus are already alternating in sign in the first difference.

The difference table algorithm will attempt to calculate the noise variance in the data, which you can display using the function show_difftab_properties_fndf,

and the output is a report that looks like

If the difference table algorithm finds at least 4 contiguous changes of sign in a difference level it is signalled as the presence of numerical error. Each level at which this occurs is displayed so that you can decide for yourself the quality of the data. The error is the estimated variance of the input data based on the variance found in the difference data at that level.

Approximate Derivatives

A simple finite difference can be used to approximate the derivative of a sample function, f(x), if the steps size, h, is small enough. If we denote the difference operator by D, the second difference by D2, ect, then

and the derivative is approximated by

Differences can be calculated with the ndiff_fndf function,

where

• hcb = the handle to the callback function
• x = the x value at which to calculate the difference
• 1 = the order of the difference
• 0.01 = the step size, h
• 'FORWARD' = type specifier for a forward difference

The HMD function for calculating the raw finite difference approximation of a derivative, as shown above, is nderivative_fndf. For example,

where

• hcb = the handle to the callback function
• x = the x value at which to calculate the derivative
• 1 = the order of the derivative
• 0.01 = the step size, h
• 'FORWARD' = type specifier for a forward difference

The sample function would be defined, for example, with

and registered as a callback function with

The problem with a raw finite difference approximation is the contention between truncation error (accuracy) and roundoff error (or bad sample function error). In order to get an accurate answer the step size, h, must be very small. But derivatives of the 3rd order and higher begin to show numerical divergence for reasons discussed above in the section on difference tables. For that reason the raw finite difference approximation has limited practical use unless you can generate samples with very high precision.

The work-around for getting a numerical derivative approximation is to create an interpolation polynomial of sufficiently high order in the interval of interest and calculate its derivative. This is what the pderivative_fndf does. An example of this function call would be

where

• hcb = the handle to the callback function
• x = the x value at which to calculate the derivative
• 1 = the order of the derivative
• 0 = the lower limit of the interpolation interval
• 2 = the upper limit of the interpolation interval
• 10 = the number of samples to take

Introduction

The standard method for the numerical integration of some function f(x) over an interval a <= x <= b is

For example, the Trapezoid Rule is

and Simpson's Rule is

where x1 = a, x2 = a + h, x3 = a + 2h = b.

The remainder term, R(z), gives an indication of the truncation error of the integration, and is usually represented as

The numerical integration library provides methods for integrating functions, f(x), using predefined formulas such as the Trapezoid Rule and Simpson's Rule. The user may also define integration formulas, inserting the necessary weights and error term factors.

An integration must be set up with three basic steps. First, a function f(x) must be defined and declared as a callback function. Second, an integration rule must be declared, which establishes a handle to be used to access the integration object. Third, the integration is performed with specified limits of integration.

The two basic function calls for implimenting an integration are

defined_rule_intg

and

integrate

the former specifies the integration formula. The latter does the integration on the limits.

Defining a Simple Integration Rule

A simple integration rule is one which defined as a weighted sum of equally spaced function samples. There are different ways of representing the formulas for numerical integration. In the following tutorial, the simple formulas conform to the representations in Moursund and Duris, "Elementary Theory and Application of Numerical Analysis."

For the Trapezod Rule,

the weights are c0 = 1/2 and c1 = 1/2.

The remainder gives an error estimate, and is composed of three factors: a truncation error factor, errFac, the mth derivative of the sample function evaluated at an unknown location, z, in the integration interval, and the step size raised to the m+1 power. For the Trapezoid Rule, the errFac is -1/12.

The rule definition is given by the code

The first argument in the call is the number of coordinate dimensions of the problem, which is 1. The rest of the arguments supply the formula parameters. After the dimension are given the integration weights, c0 and c1, preceded by the number of weights that follow. The next two arguments are not used for simple formulas. Then the magnitude of the truncation error is given, the order of the remainder derivative, and a handle to the sample function to be integrated.

The essential ingredients of the simple formula are the weights, the error factor, and the order of the remainder derivative.

The sample function, f(x), must be defined and declared as a callback. For example, the following code gives a function definition.

Then this function must be declare for use as a callback function by the integrator,

This provides a handle to the function that is passed as a parameter to the defined_rule_intg method. In practice, these sample function declarations must come before the call to defined_rule_intg.

Optionally, another argument can be given, which is a handle to a user-provided function that gives the remainder derivative of the sample function. This is needed for the computation of an error estimate. If a remainder derivative is not provided an attempt will be made to calculate a numerical approximation by finite differences.

In the special case of the Trapezoid Rule a predefined method is also provided,

which is used for simple testing.

Finally, the integration is performed. For example, with integration over the 0 <= x <= 1,

The integration formula given in mathematical tables will be in terms of multiple function samples and a step size, h, not just the two endpoints shown above:

The use of multiple function samples is equivalent to subdividing the integration interval into smaller intervals, h, and summing the results of each of the smaller integrations. The integrate method does the subdivision for you. You simply give the number of smaller intervals as the argument that follows the limits of integration. The step size is automatically calculated. For example, for 4 subintervals the function call looks like

The error bounds on the calculated integral are automatically printed to the log file. But they can be acquired programmatically by calling get_error_intg, as in the following example.

where the handle to the rule definition is again used. The value returned by get_error_intg is an array of two elements: the lower estimate and the upper estimate. These are calculated by evaluating the remainder at x = a and x = b.

There is also a predefined Simpson's Rule formula,

The rule definition for Simpson's Rule could be done with the defined_rule_intg method as follows,

Defining a Gaussian Integration Rule

A Gaussian integration rule makes use of the of the sample point locations as parameters in addition to the weights. Therefore, there are twice as many coefficients in the integration and twice the order of approximation.

An easy example of a Gaussian integration is taken from Hamming, ch 19.2.

and let

w1 = 5/9, w2 = 8/9, w3 = 5/9
x1 = -sqrt(3/5), x2 = 0, x3 = sqrt(3/5)

For the error term, see Hamming, ch 19.5.

For N = 3 (three weights),

then

The Gaussian integration rule is then defined with the defined_rule_intg command like

Here was assume that hfunc is the handle to a callback function, f(x), and hf6 is the callback handle to a literally defined 6th-order derivative of f(x). The Gaussian rule as given above may only be used on the integration limits from which the formula was derived: -1 to +1. We may therefore perform the integration on those limits like this

which is an interval from -1 to +1, not subdivided.

As long as we still use those integration limits we can also subdivide the integration interval, thereby achieving greater accuracy. Below we increase the subintervals from 1 to 4,

To be useful the integration rule must be capabale of integrating over arbitrary intervals. In order to achieve this with Gaussian rules we must call defined_rule_intg with extra parameters. After the list of x coordinates, and before the error factor, we must insert the number of integration limits and the limits themselves for which the formula was derived. For the Gaussian rule we are using for an example, the formula was derived with limits -1 to +1, so we need the three parameters 2, -1, 1 in the call to defined_rule_intg,

The Error Calculation

The default behaviour of the integrate function is to recalculate the estimate for the integration if the limits of integration are changed or the number of subintervals is changed. If you don't need the error estimate, you can turn off that automatic computation using set_intg,

where "hnd" is the handle to the integration rule. To turn it back on again, use "1" for the argument,

As was shown above, the error estimate for a simple rule is usually given by the formula

which complies with the representation from Moursund and Duris, "Elementary Theory and Application of Numerical Analysis," section 6-5.

The formula for the Gaussian error estimate, as shown above, is

and is derived in Hamming, "Numerical Methods for Scientists and Engineers", section 19.5.

This is essentially the truncated remainder from the Taylor series approximation for the integral, and is commonly called the truncation error. The mth-order derivative of the integrated function, f(x), must be evaluated at the unknown value of x = z, but we do not know the value of z (or else we would already know the exact answer). The error computaion is performed at each of the lower and upper limits of integration. So the error computation provides a good guess for the lower and upper bounds of the error correction to the result of your integration calculation.

But the error computation requires a function be given for the mth-order derivative of the integrated function, f(x). This can be defined by the user, and its callback handle passed to the defined_rule_intg function, or set later with the set_intg function. This is an optional parameter because often we do not need the error computation. But if the user does not provide an explicit function for the mth-order derivative, a finite difference approximation to it will automatically be calculated, unless you set "errflag" to 0 as discussed above.

It is often inconvenient to actually work out the mth-order derivative for a function, f(x). So the HMD integration function estimates one for you. This is done by creating an interpolation polynomial that approximates the function, f(x), within the given limits of integration. Then its mth-order derivative is evaluated. In order to obtain a reasonable degree of acuarcy for this estimate a polynomial with twice the order as the derivative is created. The default method for creating the interpolation polynomial is by construction of the Vandermonde matrix.

The most important consideration for the error calculation is that you must provide the error factor, Em, and the order, m, of the remainder derivative. For all of the well-known formulas these parameters have already been derived. For example, for the Trapezoid rule Em = -1/12 and m = 2, and for the Simpson rule Em = -1/90 and m = 4. The derivations of these parameters can be found in Moursund and Duris, Hamming, and other sources.

Introduction

ODE means Ordinary Differential Equation, more specifically a first-order, differential equation in one independent variable. A first-order ODE is usually represented, for numerical work, as

and the solution is written formally as the indefinite integral

Over a sufficiently small interval in x, x1 - x0 = h, the derivative f(x,y) is assumed to be nearly constant, and therefore the indefinite integral can be approximated by the simplest quadrature formula

Given known values of x = x0 and y = y0 at the lower limit of integration the value of the integral at a small distance h from the starting point can be approximated by the formula

which is called Euler's method. After the value y = y1 is calculated at the point x = x1, the point x1 is treated as the new starting point. The procedure is repeated to obtain y = y2 at x = x2.

Example

The ordinary differential equation y' + y = 0 is known to have the solution y = exp(-x). Using initial conditions x0 = 0, y0 = 1, and h = 0.1, the value of y = y1 is calculated:

y = y0 + f(x0, y0) * h = 1 + (-y0) * 0.1 = 1 + (-1) * 0.1 = 0.9

In the HMD interactive interpreter (hmd -i) you can type in the following statements.

Type y again without a semicolon and its value will be printed. Now press the up arrow key to repeat the last command,

You can step forward in the integration by pressing the up arrow and return keys.

The result of performing the Euler method with this equation and these initial conditions is shown in the table below, along with the accurate values in the column under Y.

x	y	Y	err = 100 (y-Y)/Y
0	1	1	0
0.1	0.9	0.9048	-0.53050
0.2	0.81	0.8187	-1.0627
0.3	0.729	0.7408	-1.5929
0.4	0.6561	0.6703	-2.1185
0.5	0.59049	0.6065	-2.6397
0.6	0.531441	0.5488	-3.1631
0.7	0.478297	0.4966	-3.6857
0.8	0.430467	0.4493	-4.1916
0.9	0.38742	0.4066	-4.7172
1.0	0.348678	0.3679	-5.2248

The HMD library contains a function for performing the Euler integration called "ode_euler".

In the function arguments, hx is the step size, x0 the initial x position, and y0 the initial y position. N is the number of steps to perform in the integration. The arguments that follow are used to specify how the ODE function y' = f(x,y) is given to ode_euler. A function f(x,y) will be defined in the interpreter workspace, and it will be used as a callback function invoked by ode_euler at each integration step. The fifth argument is the name of the previously-defined callback function. This callback function must always be defined to take 2 double-type arguments, where the first will be given the x coordinate, and the second will be give the y coordinate.

For example, the function f(x,y) as given previously, y' = -y, would be defined in the HMD workspace as

When working in an interactive session (hmd -i) it will not be possible to enter the above function definition line by line. It can be written into a file by itself, say my_fxy.hmd, and then source'd into the interactive session with the source command.

In the particular version of f(x,y) we have chosen the x variable is not used, but the ode_euler callback functions must be defined in this way.

So a complete sequence of commands in the interactive session would look like

Data

The data from the integration is returned in the "output" variable. This variable is a "structure", similar to the struct type in C. It has at least three members: data, rows, and columns. The members are accessed by following the variable name, "output", by a dot and then the member name:

output.rows
output.columns
output.data

In the example above the data member will be an array of double-precision numbers. The array will have 3 rows and 10 columns. The columns contain the integration steps, and the rows contain x, y, and y'. So, if you print the value of output.rows, it should be 3:

The values in the data array can be printed to the screen using the "dump_block" function.

Plots

The output data can be plotted using the HMD plot command. The command

will plot the x-y values in the data. You can also plot the y derivative values that were computed at each point with

The first argument to plot is the data set to be used for the x data set. Since there is no option specifying the row number it is assumed to be row 1 (the rows are counted 1, 2, 3). The second data argument is assumed to be the first y-axis data set, as there can be multiple y data sets. Because the data argument has more than one row it is necessary to specify the row number with an option string.

The HMD plot command is a front end to the Gnuplot program, which must installed already on your computer.

ODE Solvers in HMD

The HMD program provides several ODE solvers as single function calls (as opposed to building blocks, discussed later):

A trivial Euler model,

The 2nd order predictor/corrector, that is, "Heun" model,

A 4th order Adams predictor/corrector model,

The common 4th order Runga-Kutta model,

You need only supply the function y' = f(x,y) and the initialization arguments,

hx = step size
x0 = initial x value
y0 = initial y value
N = the number of steps to calculate

and for the prdictor/corrector models an additional argument, Ni, which is the number of times to iterate the corrector equation.

The return value is a block of data that can be view or saved using the dump_block function (see the Data discussion above). The block of data is a matrix with as many columns as there are steps in the calculation (and the initial values). There are rows for each order of derivative used in the model plus one for the x values. The first row are the x values. The second row contains the y values, and the third row contains the first derivative values.

Building Blocks for ODE

The ODE problem can be assembled in steps of elementary building blocks. The following functions are used for that purpose.

init_ode
set_function_ode
compute_derivatives_ode
step_ode
add_datatable_ode
get_ode_data

The elementary ODE functions are used for creating a specific ODE model not provided in the HMD library. To illustrate their use, consider the Euler model presented above.

Define the initial constants for the ODE as before.

Then "open" a handle to a new model with the init_ode function,

The init_ode function requires some extra parameters that were not introduced before. These new parameters allow ODE models to be assembled as either single-step or as multistep methods. This will be explained in the following sections. The Euler ODE method is a single-step, first-order method, and so we set the number of (past) steps to 1,

and we set the highest derivative order to 1,

The array argument at the end provides the "weights coefficients" for previous step values and derivatives. For the Euler method this is trivially 1, and 1. That means the coefficient of the last y value is 1, and the coefficient of the last y' value is 1.

The integrator must be given the name of the ODE function to call, and this is done with the set_function_ode function,

We must use, for the first argument, the handle returned from the init_ode function call to identify the model. The second argument specifies which stage of the computation will use the callback function. A "stage" is one of a sequence of equations computed for one step of the integration. For example, in a predictor-corrector model the predictor is stage 1, and the corrector is stage 2. For this example there is only one stage, so the second argument is trivially 1, meaning stage 1. The third argument specifies what derivative order to be assigned to the callback function. The Euler method has only a first-order derivative, so this is set to 1.

The set_function_ode function requires (at least) 2 parameters to specify the callback function. The last 4 arguments shown above for a standard set of callback paramaters: the callback function name, the number of arguments that will be passed to the callback, and a list of symbolic representations for the callback function parameters. $1 means x, $2 means y, $3 means y', $4 means the second derivative, etc.

The work of integration of the ODE is performed by the step_ode function. This function does all the computations in the step from position x[n], y[n], to position x[n+1], y[n+1]. The step_ode function may be called for each step, or it may be called once with a repetition parameter which will cause it to calculate N steps forward.

At the end of each step the step_ode function performs the calculation of y' = f(x,y) for the y value just computed, and its value will be used for the next step in the integration. That means that before the first step the derivative y' = f(x,y) must be manually calculated for the initial position x0, and y0. This is done with the compute_derivatives_ode function. It must be called before the first call to step_ode.

The data values from each step in the integration can be printed to the console by setting a global variable, ode_step_print, before beginning the steps.

And the output prints can be turned off again by setting it to 0.

But the step_ode function can automatically write the output data into a data table with a given number of steps. To setup a data table you must call the add_datatable_ode function before any computations are performed.

The number, N, of columns of the data table must match the number of steps to be calculated with step_ode. The number of rows of the table will be automatically assigned to match the highest-order derivative in the integration. That is, the data table will have in its first row the x values. The second row will contain the y values. The third row will contain the y' values. For higher-order calculations, such as a Taylor method, there will be more rows to contain higher derivative values.

To obtain the data in the data table as an array (matrix) in the HMD workspace, use the get_ode_data function.

The return value from get_ode_data will be an output struct with the same format as shown above for the ode_euler function, that is, it will have at least the three members

output.rows
output.columns
output.data

This data can be saved or viewed with the dump_block or plot functions.

Multistep ODE Building Blocks

A multistep ODE integration uses y values from past steps in the integration. With a choice of appropriate coefficients as factors on these past step values the order of the integration will increase. For example, the Euler integrator shown above can be rewritten as

y[n+1] = y[n] + h * y'[n]

and is a single step integration. It is single step because the indices on the right hand side select only the values from the last integration step. An integrator that uses the last two steps might look, for example, like

y[n+1] = y[n] + A * y[n-1] + h * (y'[n] + B * y'[n-1])

The integration uses the steps n and n-1 as shown on the right hand side. Thus, the second argument to the init_ode function is used. The value of np specifies the number of (past) steps to be used. The Euler model uses only one past step value, so the value of np would be 1.

Most often a multistep model will also be a multistage model, which means an integration step will be performed as two (or more) stages. One type of multistage ODE model is a predictor-corrector model. Each equation in a multistage model is constrained to be a difference equation in the past y and y' values with a choice of coefficients. So the specification of an equation to be used as a stage in the integration step will require

1. np, the number of past step values to be used in the stage
2. nO, the order of the highest derivative (usually 1) to be used
3. the coefficients of the past y and y' values

The init_ode function automatically creates the first stage in a multistage model. More stages can be added to the model using the add_branch_ode function. The equation used for a stage of the integration is also called a "branch". The add_branch_ode function requires two arguments: the handle to the ODE model that was returned from the init_ode function, and a matrix of coefficients, called "weights".

As an example, consider the simplest multistage model using only the last step in the integration. It is variously called the Heun model, or the implicit Euler model, or a Backward Difference Formula (BDF).

y[n+1] = y[n] + (h/2) * (y'[n] + y'[n+1])

The right hand side contains, implicitly, the value of the derivative at the n+1 step. But the n+1 step has not been calculated yet. It is therefore called an implicit method. A starting value for y'[n+1] must be obtained first. The stage that calculates the starting value is called the predictor, and the final, implicit, calculation stage is called the corrector.

The first stage of the integration is performed using the formula

y[n+1] = y[n] + h * y'[n]

which is just the euler calculation. And the trial value of y'[n+1] is calulated with f(x+h, y[n+1]). The init_ode function sets up the first branch

The coefficients for the first branch are given by the 2-row, 1-column matrix

1.0
1.0

Row 1 of the weights matrix gives the coefficients on the past y values, and row 2 gives the coefficients on the past y' values.

Each branch of the ODE model is given its own handle number to be referenced later. The init_ode function always assigns a handle equal to 1 to the first ODE branch.

The corrector branch of the ODE

y[n+1] = y[n] + (h/2) * (y'[n] + y'[n+1])

is added using the add_branch_ode function.

Here the coefficients assigned to this branch are given by the matrix

1.0
0.5

because the coefficient on y[n] is 1, and the coefficient on y'[n] is 1/2.

The coefficient on the implicit value, y'[n+1], is provided in the next argument, the row array

0.5

0.0

The second entry in this array, as well as the last parameter [1,0], will be explained later. For two-stage predictor/corrector models the last two arguments of add_branch_ode will always look like

...,

[B, 0]

[1, 0]

where "B" is the implicit coefficient on y'[n+1].

The size of the weights matrix must correspond to the ODE parameters given in the init_ode function call. The matrix will have np columns and nO+1 rows. The format of the matrix is shown below,

	n	n-1	n-2	...	(past steps)
y	A1	A2	A3	...	(0th derivative weights)
y'	B1	B2	B3	...	(1st derivative weights)
y''	C1	C2	C3	...	(2nd derivative weights)
.	.	.	.	...

Note that each branch, or "stage", of the model must be provided with a different weight table as shown above. The weights for the first branch are always provided as the last argument to the init_ode function, and for each successive branch in the second argument of add_branch_ode.

The following example shows the code for building the Heun, or "implicit Euler" model.

Notice the extra function

in the above example. By setting the 'bwmcnt' parameter to 3 the y'[n+1] value will be reiterated 3 times to improve the accuracy.

There are some details you should observe. Notice that the size of the data table created with add_datatable_ode is equal to the number of steps given in step_ode. You must include the call to compute_derivatives_ode before the call to step_ode so that the table of past steps will be initialized. The call to set_ode with 'bwmcnt' is not necessary, but it is known that reiterating the corrector will slightly reduce the truncation error of the step. The accuracy does not improve after about 3 or 4 iterations. Lastly, observe that the single call to step_ode could be separated into repeated single steps with the call

Runga-Kutta Building Blocks

Runga-Kutta is the name for an approach to solving ODEs. It is usually called a "single step" approach, as opposed to the multistep approach of predictor/corrector models. That means that a Runga-Kutta model does not need to include the results from past step calculations to calculate the next step. The most well-known Runga-Kutta calculation is the 4th order method, given by the following algorithm:

y(n+1) = y(n) + (1/6) * (k1 + 2*k2 + 2*k3 + k4)

k1 = h * f(xn, y(n)) k2 = h * f(xn + h/2, y(n) + k1/2) k3 = h * f(xn + h/2, y(n) + k2/2) k4 = h * f(xn + h, y(n) + k3)

It is not a multistep calculation, but it does use multiple stages. These stages (k1, k2, k3, and k4) are obviously of a different nature than the two stages of predictor corrector models. But there is a similarity from the computer programming point of view (see The General Linear approach, below).

The building block approach to Runga-Kutta models follows the same pattern of assembly as for the predictor/corrector models. A "branch" of the ODE will be created for each of the calculation stages k1, k2, k3, k4, but with a minor difference that a special fifth stage is supplied. Also, because the Runga-Kutta models calculate a weighted sum of the k's, an extra building block function called a "quadrature" must be provided.

An example set of building blocks that computes the 4th order Runga-Kutta models is shown here.

The correspondence of the k stage calculations to the weight arrays given as arguments to the branches will be preesented here. But a more complete picture of this algorithm must wait for the discussion on The General Linear ODE model below.

The first branch of the ODE is created by calling the init_ode function, but for a Runga-Kutta model the first branch does not really start the calculation. It just provides y0. The calculation of the k's starts in the second branch and finishes with the special fifth branch. Think of it this way: whereas in the predictor/corrector model the result of a branch is the value y[n+1], a branch of the Runga-Kutta model produces the "y" input to f(x,y) used in the next k calculation. Thus we have

y(n) output from branch 1 y(n) + k1 / 2 output from branch 2 y(n) + k2 / 2 output from branch 3 y(n) + k3 output from branch 4

and k4 itself needs the special fifth branch to be calculated. The calculation of each k(n)/h is performed as the f(x,y) using the output of the previous branch for the y parameter. The quirk of this algorithm is that the k(n) values are logically matched according to the "stage" of the calculation, i.e., k1 to stage 1, k2 to stage 2, etc. But k1 is actually computed at the beginning of branch 2, k2 at the beginning of branch 3, etc.

The arrays [1; 0] supplied as the second argument to add_branch_ode are coefficients of 1 to y(n) and 0 to y'(n). This is what produces the y(n) part of the branch outputs shown above.

The third argument to add_branch_ode are the "implicity" weights with the following correspondence

branch 2: [1/2, 0, 0, 0] --> coefficient of k1 branch 3: [0, 1/2, 0, 0] --> coefficient of k2 branch 4: [0, 0, 1, 0] --> coefficient of k3

The last argument to add_branch_ode gives the weights that produce a fractional step size in the variable x as the first argument to f(x,y).

The inclusion of add_quadrature_ode after all the branches have been created is what adds the final calculation of the Runga-Kutta step,

y(n+1) = y(n) + (1/6) * (k1 + 2*k2 + 2*k3 + k4)

where obviously the weight array gives the coefficients of the k(n) in the weighted sum.

The General Linear ODE Approach

John C. Butcher in his book, "Numerical Methods for Ordinary Differential Equations", Second Edition, presents the General Linear Method for representing ordinary differential equations. The most significant result of the General Linear Form is to show a correspondence between the computations of Predictor/Corrector models and Runga-Kutta models.

The General Linear Form is a matrix M given as

The submatrix U contains the coefficients of past step values, the submatrix A contains coefficients of h*f(x,y) that will be used to calculate either y'(n+1) for a predictor/corrector, or the k's for a runga-kutta model. The submatrices B and V contain coefficients for the final, combining sums that result in y(n+1).

The solution process is represented as a matrix equation

G = M * X

or

This can also be written

The size of the input vector is determined by the "s" stages of the calculation and the number "r" of past step values used. This complicated- looking set of equations can be made easier to see by comparison to the known ODE methods. The first set of equations, Y_i, represents the two equations in the predictor/corrector method. Or it represents the calculation of k1, k2, k3, k4 in the 4th order Runga-Kutta method. The second set of equations, W_i, represents the weighted sum that finishes the Runga-Kutta step.

To see this better, consider the 2nd order predictor/corrector step calculation,

The first calculation contains only the second sum of Y_i, that is, all the coefficients a_ij are zero. The second calculation of the predictor/corrector has the input, F_1, as the f(x,y) value using as y the Y_1 from the previous stage. This could also be shown as

In this example the submatrix A would be

and the submatrix U would be

In this case, the single equation for W_1 is just a repetition of the Y_2 equation. This is just saying that the predictor/corrector needs no final weighted sum. But formally, to satisfy the matrix equation we could write

From the above example it should be obvious that the General Linear procedure is not as simple as a linear algebra calculation. The input vector must be updated after each row computation, as can be seen with the F_i values. But there is another piece of information missing in the representation: the step position in x where f(x,y) is to be evaluated for the value of F_i.

As another example, consider the 4th order Runga-Kutta step calculation,

This would be represented in terms of the General Linear matrix as

In terms of the linear equations for Y_i and W_i this would be written as

It will be seen that k1 = h * F_1, k2 = h * F_2, etc. As mentioned above, the F_i are updated after each row calculation. So the k1, k2, k3, k4 are calculated between stages. The k4 is calculated after stage 4. Also, observe that Y_1 is the input to F_1, Y_2 is the input to F_2, and so on, as required for the calculation of the k values.

In the W_i calculation the F_i are NOT updated after each row. Now the F_i are all the values computed for the Y_i equations.

The HMD Implimentation of the General Linear Form

The HMD implimentation of the General Linear ODE approach obeys the following rules,

1. Each row of the upper partition of the General Linear matrix is represented by an ODE "branch" and is roughly equivalent to a "stage" of a multistage model. So the Y_1 calculation is set up with the init_ode function, the Y_2 calculation is set up with the add_branch_ode function for the second branch, and so on.


2. Each row of the lower partition is represented by a "quadrature" calculation that is set up by calling the add_quadrature_ode function.


3. The result of the last quadrature calculation will be the output for the ODE step. If no quadrature calculation is defined then the output of the ODE step will be the result of the last branch calculation.


4. The coefficients in the U submatrix correspond to the "weight table" that is given as a parameter to each of the add_branch_ode function calls. Each row of U corresponds to one branch weight table as a concatenation of its columns. For example, if the weight table, w_ij, has 2 rows and 3 columns then that row of U would correspond to w_11, w_21, w_21, w_22, w_31, w_32. The coefficients of the A submatrix are given by the "implicit weights" that are passed to add_branch_ode. The weights arrays (matrices) given on each call to *add_branch_ode* constitute a single row of the A or U submatrix of the General Linear matrix.


5. The 4th parameter to the add_branch_ode function call is an array that has no correspondence with the General Linear matrix. It provides the fractional positions within an ODE step where the F_i = f(x,y) will be calculated. For a predictor/corrector model this is trivially 0 or 1 (x0 or x0 + h). For a Runga-Kutta model it is 0, 1/2 and 1 (x0, x0 + h/2, x0 + h).

The core of the HMD integrator for ODE is the internal Ode_Step function, which first loops through the ODE branch objects (stages) then loops through the ODE quadrature objects (weighted sums).

Stopping ODE Integration with a Trap

The need might arise where the integration should be stopped if the integration state, y, takes values that are out of bounds with respect to the derivative calculation. For example, the value of y becomes less than zero or becomes greater than some upper bound. For this circumstance you can set up a "trap" with the function call trap_ode_model.

The function syntax requires four parameters: the model handle,an integer specifier of the order of derivative of the ODE state to trap, a string specifier for the type of test condition, and the value that it should be compared to. The integer specifier for the order of derivative starts at 0 for y, 1 means y', 2 means y'', etc. The string that specifies the comparison type is a two letter, upper case identifier similar to the Fortran test condition:

• LT, less than
• LTE, less than or equal
• EQ, equal
• GT, greater than
• GTE, greater than or equal

Set up a System of Equations

Systems of ordinary differential equations are integrated by setting up each equation as an individual model, as above. Then some extra functions are called to "glue" the separate models into one system.

You initialize the system of equations by providing a list of the model handles that will become part of the system. This is done with the init_ode_system call, which returns a handle to the new system of equations. An example with two equations would look like

With a system of equations there must be an extra configuration of the callback function arguments for those functions that define the first derivatives. This is already done for each model setup with the call to set_function_ode, so it seems reduntant to have to do it again when combining the equations into a system. The reason for this is that by default for a single ODE the list of function arguments in set_function_ode is assumed to define x, then y, then y', y'', etc. That is, the argument placeholders of $3, $4, and above are considered to specify higher derivatives as arguments. But when defining a system of equations some of these arguments must be state variables for the cointegrated equations. For example, for a system of two equations, with states y1 and y2, the argument list for any one of the derivative functions could be x, y1, y2.

The second argument list specification for each of the ODE models is performed with the set_function_args_ode_system call. An example of using this function for a system of two equations in which the first model (equation) takes x, y1, and y2 as parameters would be

Notice that the last parameter to set_function_args_ode_system is a string in single quotes containing $, $1_0, and $2_0. The $ symbols are use in a slightly different format than in set_function_ode. The state variables are numbered 1, 2, 3, ... for y1, y2, y3, etc. The number following the underscore specifies the order of the derivative. So $ means x, $1_0 means y1, and $2_0 means y2.

The first 4 arguments to set_function_args_ode_system are

1. ODE system handle
2. ODE model handle
3. ODE branch handle (usually just "1")
4. The order of the derivative function (usually just "1")

For a system of ODE you do not call step_ode on individual ODE models (individual equations). Rather, you must call step_ode_system.

After the integration is complete you may retreive the data from the ODE models individually, as if they were single ODE models. You can also get the data in the form of a system output by calling get_ode_system_data. For example,

In the above example the variable outstates will be a struct with the same format as an output struct for a single ODE model. It will have a rows member, a columns member, and a data member. The first row of the data member will be the x data; the next row will be the y1 data, the next the y2 data, and so on.

Starting a System of Equations with Another

A multistep, predictor/corrector model that is greater than second-order needs to be started with another lower-order predictor/corrector model or a Runga-Kutta model.

It is a simple matter to join the two models so that the main model starts where the starter model stops. The two models must share one "step table." Each model has a step table, which is an internal matrix of past ODE step calculations and their derivative calculations. The step table has as many columns as needed to store the past values. It will have as many rows as are needed to save the highest order of derivative needed by the model (these are specified in the init_ode call). The starter model needs a smaller step table than the main model. So the main model keeps its step table, and the starter model will be set to use the main model's step table.

First perform the setup of both models, but do not call the step_ode integrator yet. Then call share_step_table_ode.

The first argument, h_model1, in the above code example is the handle to the starter model. The third argument, h_model2, is the handle to the main model. Arguments two and four should be zero (These are actually meant to be handles for specific ODE branches within each model, but that usage would entail a somewhat complicated integrator).

When integrating a system of equations we must call share_step_table_ode for each model (equation) in the system. For example, if the system is a set of three equations we must call share_step_table_ode three times. Model 1 in system 1 will be matched with model 1 in system2; model 2 in system 1 will be matched with model 2 in system 2, etc.

But a system of equations, as opposed to a single model, will require an additional function call for gluing the starter system with the main system. We must call share_states_table_ode.

In the above function call the first argument specifies the starter system, and the second argument specifies the main system.

Introduction

The finite element method is an engineering method for solving two kinds of equations

where the first equation represents static forces on a system, and the second equation represents vibrational modes. In a nutshell, the first equation is the starting point for defining static problems, and the second for defining dynamic problems. The mathematics for diverse areas of engineering, such as fluid mechanics, can be shoe-horned into forms that resemble the above equations.

The solution of these partial differential equations is a quadrature (definite integral) as opposed to an iteration (indefinite integral) as is performed on ODE problems. This is possible because the finite element method is always applied to boundary-value problems as opposed to initial-value problems. You approximate the funtion phi over a small region as an interpolation function in two or three dimensions. Rather than working with polynomial coefficients the interpolation is transformed into node values and influence functions. Then apply the Laplacian operation (two differentiations), and perform the quadrature.

The result is a numerical solution over a tiny region, i.e., an element. For a single element the node values and the boundary values are equivalent. In other words, nothing is solved. The usefulness of the method is achieved by combining many elements into a set of equations. The outer nodes are known boundary values, and the inner nodes are simultaneous solutions of the problem.

This simultaneous solution results in a stiffness matrix, S, and the forces end up in a vector, f.

The simultaneous solution of the second equation results in an inertia matrix, M, and the second equation solves for the mode frequencies and the mode shapes of a vibrational problem.

FE Model Components

The finite element model is constructed from a mesh file and an HMD script file.

The mesh file contains the partition of some geometric region into elements that are defined by vertices and their locations. Nodes are assigned to the vertices and are shared by neighboring elements. The mesh file is created with the GMSH program, which is freely available from the GMSH web site. Documentation on how to use GMSH is found there.

The HMD script reads the mesh file and creates internal data structures for the nodes and elements. During the creation of the mesh file with GMSH every element is designated as belonging to a physical region and given a region number. In the HMD script the region numbers are assigned material coefficients, boundary values, and forces (that is, applied influences). The HMD program assembles the matrix equations from the nodes and elements and solves the matrix equations.

An output file is generated for each finite element model. Although the output data can be viewed and manipulated with any program you like, it is expedient to use the GMSH program for working with the output. For that you need the HMD data file and the GMSH mesh file. There are HMD functions that can produce data files in the GMSH post processing format. You can then use the GMSH display to examine your model.

The Simple FE Function

HMD provides a simple, one-function call approach to the static finite element equation defined above.

where

output = a struct with (at least) three members: rows, columns, data. The data member will be a block of double-precision numbers giving the nodal values.

keyword = a two-letter text string: 'bc', 'gr', 'ms', 'fc'.

• 'bc' means boundary condition value
• 'gr' means the stiffness value
• 'ms' means the inertia or mass value
• 'fc' means the value of the forcing term.

region_number = an integer identifier of the region to be given a value. Regions are created during the creation of the mesh file.

region_value = a double-precision number to be assigned to the region.

Example,

Given that you have created a mesh file called triGreen.msh, call

where we used the dump_block function to write the output data to the file triGreen.dat.

The output from an HMD finite element model can be displayed using GMSH. An AWK script provided in the HMD share directory (/usr/local/share/hmd/scripts) called gmsh2post.awk which can create a GMSH post-processing file. For a model based on the following mesh

The output displayed with GMSH looks like

Interface with GMSH

Documentation and tutorials for using GMSH are available at the GMSH web site. The user must first create a geometry file, building points, lines, surfaces, and volumes. It has a graphical user interface for building the geometry visually. Or the geometry file can be created from a text file. The following text fiel was used to create the geometry for the above model.

If you use the graphical interface the result will still be a text file as above. Note, however, the comment in the file

gmsh -2 -format msh1 tri_green.geo

The mesh must be created with format "msh1", that is, with the legacy GMSH mesh format. The HMD software development began when GMSH was still a relatively new program. All HMD mesh interfaces use the legacy GMSH formats.

The HMD program can send native GMSH commands to the GMSH program over a socket connection. The GMSH program must be started already in interactive mode using the "listen" option:

When running HMD you must tell it what socket name to use when establishing a connection to GMSH. This is done with the HMD -g option. The default socket name used by GMSH will be .gmshsock located in the home directory for your user account. So HMD would be involked as

The simplest way to send commands to GMSH is with the

function call. If a socket connection is not already established an attempt will be made to make the connection before sending the GMSH command. Then the socket is left open for the remainder of the HMD session.

If the socket connection to GMSH is not explicitly closed during the HMD session the HMD program will attempt to close it when it exits.

The GMSH socket can be opened manually in the HMD script with the function

Or, when opening and closing the socket multiple times the socket name argument may be omitted on subsequent open calls,

in which the last socket name used will be used as the socket name. Or, it the first open call, and a socket name was given on the command line, that socket name will be used.

Closing the socket connection is performed using the function call

Note that no argument is needed--no socket handle is needed. That is because there can be only one socket connection to GMSH open at any one time. For example, if you try to call gmsh_open with a different socket name when a socket session is already open, the open call will be ignored until the current session is closed.

Low-Level FE Commands

For creating more specialized finite element models than can be performed with the simple FE function, there are low-level finite element commands available. Documentation of these commands is a work in progress. The following steps outline the creation of a static finite element model:

1. Declare the mesh

   > mesh smesh;



2. Set the type of mesh

   > smesh.meshtype = GMSH;

   
   Notice that the name of the mesh is now treated like a struct or class object.


3. Set the filename of the mesh to read

   > smesh.filename = 'tri_green.msh';



4. Declare the model.

   > model tri, 'smesh';

   
   The first argument is the name of the model. The second argument is the name of the mesh that was already declared.


5. Set the type of model

   > tri.modeltype = FINITE_ELEMENT_E;

   
   The model types are FINITE_ELEMENT and FINITE_ELEMENT_E. The difference is that FINITE_ELEMENT is the old algorithm.


6. Set the equation type

   > tri.equation = POISSON;



7. Read the mesh file and load into memory

   > loadmesh smesh;



8. Set up special model flags

   > createsparsity 'tri', 0x00;



9. Set the boundary conditions. Use the "setregionbc" and "setregionbcn" commands for any region that should be set. Arguments are model name, region number, and boundary value.

   > setregionbc 'tri', 1, 0.0;

   
   which sets the boundary value of region 1 to 0.0. The difference between commands is that "setregionbc" is for Dirichtlet and "setregionbcn" is for Neumann BC's.


10. Set the stiffness coefficients in any region. Arguments are model name, region number, and stiffness value.

    > setregiongrad 'tri', 2, 1;

    
    which sets the stiffness value in region 2 to 1.0.


11. Set the "forcing", or "driving" coefficients. For any region required use the command "setregionforce". Arguments are model name, region number, and forcing value.

    > setregionforce 'tri', 2, 1;

    
    which sets the forcing cefficient of 1 in region 2 of model named 'tri'.


12. Assemble the FE matrix

    > buildglobalgrad 'tri';



13. Assemble the vector (the "right-hand side")

    > buildglobalvector 'tri';



14. Assign the matrix just assembled to be the stiffness matrix for the model,

    > setdivergence 'tri', 'tri-grad';

    
    The second argument is the default name given to the matrix when it was created by "buildglobalgrad".


15. Assign the vector just created to be the forcing vector for the model,

    > setsource 'tri', 'tri-force';

    
    The second argument is the default name given to the vector when it was created by "buildglobalvector".


16. Insert the boundary values into the matrix and the vector,

    > insertbcmodel tri;



17. Solve the matrix equation--that is, solve for the scalar values at each of the nodes,

    > solve_div 'tri', 'tri_green.dat', 0x00, 'all';

    
    which will write the output to the file 'tri_green.dat'.


18. Use the EXIT command to finish the script.

Introduction to Gnumbers

A "Gnumber" is a number type in HMD for performing integer-only, rational number arithmetic of arbitrary-length precision. A Gnumber looks like

123 / 7

whereas a normal, standard HMD number is a floating-point number and looks like

17.5714

A Gnumber can represent a whole number of arbitrary length. A machine number will always be limited in length. For instance, a 32-bit integer can have at most 10 digits. A Gnumber integer can have any number of digits up to the limits of available computer memory. But to create an integer of more than machine digit length you must enclose it in single quotes.

You have to create a Gnumber in a special way, that is, you cannot just enter an integer or divide two whole numbers like you can in Lisp or Maxima. You execute the gnumber function call which creates a variable containing a Gnumber.

A calculation performed with Gnumbers will never be reduced to a floating-point number, so in some sense the result of a calculation remains "exact." Divisions result in a rational number, composed of integers, reduced to their lowest factors.

Computation with Gnumbers is many times slower than with ordinary computer types such as double-precision floating point numbers. Gnumbers are software objects and fundamental arithmetic is performed in software, whereas a floating-point number is a machine object, and computations are performed by the processor. Therefore, Gnumbers are not of much interest for numerical modeling. Gnumbers would of interest in number theory, or for testing the limits of truncation error for a certain type of calculation, or other special purposes.

Another feature of a Gnumber is that it can be defined with any desired radix up to a maximum value. But computatiuons between any two Gnumbers must be performed with the same radix. Internally, a Gnumber is a list of numerical digits, where each digit is a natural number between 0 and radix-1. For example, in base 10 a digit can be between 0 and 9. Each higher-order digit represents another higher power in the radix. A radix can be chosen as any natural number between 2 and 32768, with the radix 0 having a special significance. When a Gnumber is printed any digit with value greater than 9 will be represented by more than one decimal digit and separated from others by vertical bars.

Defining a Gnumber

A Gnumber is created with a function call (like a constructor) and lives in the form of a software variable only. There is no literal form of a Gnumber. To create the number 2 as a base 10 variable you would type

or to create the number 2 in base 7 you would type

Although it will be printed in human-readable form, there is no literal representation for creating a Gnumber. In some other computer languages, say Lisp or Maxima, you could just enter numbers without decimal points, like 10 and 7, and get the rational number 10/7. Such computer langauges decide whether the number should be integer or floating point by whether it finds a fractional part after the decimal point. In HMD all literal numbers are basically floating-point double-precision numbers whether or not there is a decimal point. To get a pure rational number you must define it with the gnumber() function call. So, the number 3/2 in base 10 is created with

After the radix, each position in the gnumber() function call places its integer into a specific location of the resulting Gnumber:

(1)			(3)
-----	+	i	-----
(2)			(4)

So, creating a Gnumber with

results in

There may be a total of five arguments to the gnumber function call:

gnumber( radix, numerator, denominator, imaginary numerator, imaginary denominator );

The first two arguments are required, the last three are optional.

Rationalization

The act of rationalizing a floating-point number mean converting it to a ratio of whole numbers. For example,

Complex Integer Numbers

A Gnumber can have a "complex" part. Of course, this is a stretch of the definition of complex numbers. True complex numbers are a generalization of the set of real numbers, and Gnumbers are only rational numbers. But an approximation of complex numbers using rational numbers only can be created by entering numbers into the last places of the gnumber function:

Also, a floating-point version of a complex number can be entered into the gnumber function to create a complex Gnumber:

Arithmetic with Gnumbers

There are 5 basic operations defined for Gnumber symbols: addition, subtraction, multiplication, division and exponentiation. These are performed with the usual arithmetic operators +, -, *, /, and ^. For example,

A Gnumber may be raised to an integer power (possibly negative) by another Gnumber. For example,

The exponent operator will raise an error if you try to raise a Gnumber to a fractional power. The reason for this is that a non-integer exponent on a number is a hybrid concept that partakes of both mulltiplication and root finding. You can raise x to the p/q power by first raising x to the p power, then taking the qth root of the result, or in the opposite order. The result in either case is usually the same, but the decision is yours, the mathematician. The nth root of a Gnumber can be calculated with the nth_root_number function call,

which uses a recursive Newton-Raphson algorithm.

Change of Radix

The first argument to the gnumber creation function is the desired radix, or base of representation. If you have made computations that result in one Gnumber of a given radix, you might want to see the same number converted to a different radix. This can easily be done by generating another Gnumber using the first as the input argument. Beginning with one Gnumber, say, g1, it would be converted to g2 with

The use of an arbitrary radix greater than 9, while of academic interest only, is included merely because it was an unexpected outcome of the Gnumber algorithm design. Gnumber arithmetic operations work digit by digit, like you would work it out by hand, independent of the value of the radix. The only consideration that makes a difference is how to represent it when it is printed to the screen. In a radix greater than 9 some "digits" require 2 or more placeholders (hexidecimal gets around this by extending the symbol set with A THROUGH F). In order to insure that the digits are properly separated HMD will print multi-symbol "digits" surrounded by vertical bars (the "|" symbol). For example, the number 37812 in base 60 would be created and printed as shown below,

A special radix value of 0 signifies the maximum radix, 32768. Any Gnumber having this maximum radix will automatically be converted to a base 10 number before being printed.

Using a Gnumber where each "digit" represents a higher power of 32768 is a more efficient way to carry out calculations than those using, say, base 10. The process of doing arithmetic by stepping, recursively, through each of the digits of two numbers requires software overhead hundreds or thousands of time greater than would be needed by the computer processor. A huge radix like 32768 allows many, if not most, numbers of interest to be represented by only a few "digits" of the Gnumber. A computation that operates on only one or two digits will be much faster than those operating on many. This is basically the internal design for how HMD accomplishes its arbitrary-precision integer arithmetic.

Limitations

Entry of numbers into the gnumber() function are limited in size to that allowed by the IEEE double-precision number size. In order to get a huge number at gnumber creation time you will have to enclose the number in single quotes. For example, entering 142857142857 will result in a malformed number (which will probably show up negative). But entering '142857142857' will result in a valid Gnumber. Any number enclosed in quotes is not allowed to have a fractional part or an imaginary part.