Tutorial Introduction to Aladdin

RUNNING ALADDIN

ALADDIN's command line arguments are setup in a very flexible way.

The program can be started by typing ALADDIN, the name of the executable program, followed by zero or more command line arguments. A list of ALADDIN command options can be obtained by typing either

    ALADDIN 

by itself, or

    ALADDIN  -h

MODES OF INPUT

Two modes of operation are supported:

Input from the Keyboard

    ALADDIN -[ks]

Input from a File

    ALADDIN -[sf]   

Command options have the following meaning:

The command ALADDIN must be accompanied by either the -f flag, or the -k flag. The -s flag is optional. Command line options can be typed in arbitrary sequences.

Recommendation : For anything other than the most basic problem, we recommend that you prepare your problem descriptions in an input file. For example, the command

    ALADDIN -f input-highway-bridge >! output-highway-bridge

will instruct ALADDIN to read blocks of statements from a file called input-highway-bridge and redirect the program output to a file called output-highway-bridge.

Use of an input file will allow for the incremental solution to a problem.


LEAVING ALADDIN

The command to exit from ALADDIN is

   quit;

Don't forget the semicolon.

COMMAND LANGUAGE FORMAT

The ALADDIN command language corresponds to individual statements, and sequences of statements.

     statement 1;    
     statement 2;    
     ........
     statement N;

Each statement ends with a semi-colon (;).

Comment Statements

As with C, comment statements are enclosed between /* .... */ , as in

    /*
     *  =========================================================
     *  Here is a comment that introduces a block of N statements
     *  =========================================================
     */

     statement 1;         /* the  first ALADDIN statement */
     statement 2;         /* the second ALADDIN statement */
     ........
     statement N;         /* the   n^th ALADDIN statement */

Last Statement : The last statement in an ALADDIN session should be

    quit;

DATA TYPES

ALADDIN support three data types, character strings, physical quantities, and matrices of physical quantities. Here are some examples of each type:

Character Strings

A character string is simply a sequence of characters enclosed between quotes. To print the character string

     "My First Character String"
just type something like
     print "My First Character String\n";
The end-of-line character (\n) forces output onto a new line.

Physical Quantities

A physical quantity is a number plus a set of physical units. For example, the script

     2 cm/sec;
defines a velocity of two centimeters-per-second. In ALADDIN, numbers are simply physical quantities without dimensions.

Matrices

In ALADDIN a matrix is a rectangular array of physical quantities. For example, the statement

    A = [ 1, 2, 3;
          4, 5, 6 ];
defines a (2x3) array of numbers. The script
Spectra = [ 0.0 sec,  10.0 cm/sec/sec;
            0.1 sec,  12.3 cm/sec/sec;
            0.2 sec,  12.3 cm/sec/sec;
            0.4 sec,   5.0 cm/sec/sec ];
defines a (4x2) array of physical quantities. Entries in the first and second columns have units of time and acceleration, respectively. The command
PrintMatrix ( Spectra );
generates the output
MATRIX : "Spectra"

row/col                  1            2   
        units          sec      m/sec^2   
   1            0.00000e+00  1.00000e-01
   2            1.00000e-01  1.23000e-01
   3            2.00000e-01  1.23000e-01
   4            4.00000e-01  5.00000e-02

KEYWORDS

Version 1.0 of ALADDIN has only 10 keywords, namely:
     break,        else,         for,      if,    print, 
      quit,  SetUnitsOn, SetUnitsOff,    then,    while

SI and US Units

The following names are reserved for SI and US units:

SI Units   micron, mm,	cm, dm, m, km, g, kg, mg, N, kN, kgf,
           Pa, kPa, MPa, GPa, deg_C, DEG_C, Jou, kJ, Watt, kW,
           Hz, rpm, cps 

US Units   in, ft, yard, mile, mil, micro_in, lb, klb, ton,
           grain, lbf, kips,  psi, ksi, deg_F, DEG_F, gallon,
           barrel, sec,  ms,  min,  hr, deg, rad 

Some Remarks


BUILTIN CONSTANTS

ALADDIN's built-in constants are:

    Name                Value                            Purpose
    ============================================================

      DEG      57.29577951308      Radians to degrees conversion
        E      2.718281828459
    GAMMA      0.577215664901
       PI      3.141592653589

Example 1 : The script of code

    print "1 radian = ", DEG, "degrees \n";
    quit;
generates the output
    1 radian =       57.3 degrees

Example 2 : To compute the area of a circle, we coupld write:

    radius = 3 cm;
    area   = PI*radius*radius;

    print "circle radius = ", radius, "\n";
    print "circle area   = ",   area (cm^2), "\n";
    quit;
The program output is:
    circle radius is =          3 cm 
    circle area is   =      28.27 cm^2 

PRINTING STRINGS

Basic output of character strings and physical quantities is handled with the print command.

Example 1 :

    print "Here is one line of output \n";

gives

    Here is one line of output

Character strings are enclosed within quotes (i.e. "....." ).

The escape character "\n" forces output onto a new line.

Example 2 : A single print command can generate multiple lines of output. The statement:

    print "Here is line 1\n" , "Here is line 2\n";

gives

    Here is line 1
    Here is line 2

Example 3 : The tab escape character, \t, can be used to align print output into tidy columns. The single statement

    print "One \tTwo \tThree \n" , "Four \tFive \tSix\n";

will generate something like:

    One    Two    Three
    Four   Five   Six

ALADDIN makes extensive use of units as a means of clarifying matrix and finite element problem descriptions, and program output.


CONSTRUCTION OF ENGINEERING QUANTITIES

In ALADDIN, an engineering quantity is simply a numerical value followed by a set of units. The syntax for defining an engineering quantity is

    number units

Here number is an integer or floating point number (e.g. 3, 3.0 or -3.1415926), and units

    units = [Length]^a . [Mass]^b . [Time]^c . [Temp]^d

is a product of primary length, mass, time, and temperature units. The coefficients a, b, c, and d are exponents.

The multiply, divide, and exponentiation operators (i.e. *, * and ^) are used to build sets of units that are combinations of primary units.

Example 1 : The three statements

    100 m^2;            2 m/sec;            2 kg*m^2;

define from left-to-right, an area 100 meters squared, a velocity of 2 meters per second, and a mass moment of inertia of 2 kilograms meters squared.


SYSTEMS OF UNITS

Primary units may be specified in both the SI and US systems.

SI System of Units : The primary units for the SI system are:

   SYSTEM     BASIC UNIT                  SYMBOL           EXAMPLE
   =====================================================================

   SI         Length : meter              m                2.3 m;

                       kilometer          km               0.023 km;
                       centimeter         cm               100 cm;
                       millimeter         mm               1500 mm;

              Mass :   kilogram           kg               1000 kg;

                       gram               g                1500 g;
                       megagram           Mg               1 Mg;

              Time  :  second             sec              10.0 sec;

                       millisecond        ms               500 ms;
                       minute             min              60 min;
                       hour               hr               1   hr;

              Temperature : centigrade    deg_C            100 deg_C;

US System of Units : The primary units for the US system are:

   SYSTEM     BASIC UNIT                  SYMBOL           EXAMPLE
   =====================================================================

   US         Length : inch               in               12 in;

                       foot               ft               1 ft;
                       yard               yard             80 yard;
                       mile               mile             1 mile;

              Mass   : pound              lb               2240 lb;

              Time   : second             sec              10.0 sec;

                       millisecond        ms               500 ms;
                       minute             min              60 min;
                       hour               hr               1 hr;

              Temperature : fahrenheit    deg_F

Supplementary Systems of Units : Planar angles measured in radians and degrees are supplementary units, and ALADDIN 2.0 will deal with them in a consistent way.

   SYSTEM        SUPPLEMENTARY UNIT       SYMBOL           EXAMPLE
   =====================================================================

   US and SI  Planar Angle : radian       rad              2*PI rad;
                             degree       deg              360.0 deg;

Note : ALADDIN stores engineering quantities as unscaled values with respect to a set of reference units. By default, all quantities are stored internally in the SI system of units.


QUANTITY CONSTANTS AND VARIABLES

A quantity constant is simply a number followed by a set of units.

Example 2 : The statement

    2 m 

defines the quantity constant two meters.

A quantity variable is simply a quantity constant assigned to a variable name.

Example 3 : The statement

    xAccel = 2 m/sec/sec;

assigns acceleration 2 m/sec^2 to the variable xAccel .

Note : Like most programming languages, the ALADDIN language has keywords and constants whose names are reserved for a special purposes -- keyword and constant names should not be used for variable names.


PRINTING QUANTITIES

Basic output of quantity constants and quantity variables is handled by the print function.

Example 4 : The statement

    print "The x coordinate is ", 2 m, "\n";

generates the output

    The x coordinate is          2 m

In Version 1.0 of ALADDIN, the numerical component of a physical quantity will be written in the format 10.4g (this is a C programming convention). Roughly speaking, the g conversion specification will print a quantity as an integer when at all possible. Otherwise, the quantity will be printed as a floating point number, and if needed, in exponential format. Here are some examples:

Example 5 : The statement

    print "Acceleration due to gravity is ", 32.2 ft/sec/sec , "\n";

generates the output

    Acceleration due to gravity is         32.2 ft/sec^2

CASTING QUANTITY OUTPUT

The command option

           quantity ( < units > ) 

enables the printing of quantities with a desired scaling of units.

Example 6 : The script of code:

    x = 30.5 m;
    print "distance = ", x (ft), "\n";

generates the output:

    distance =      100.1 ft

Example 7 : The script of code:

    gravity = 9.81 m/sec^2;
    print "gravity = ", gravity (ft/sec/sec), "\n";

generates the output:

    gravity =      32.19 ft/sec/sec

TEMPERATURE

Temperatures may be stored and manipulated in the SI and US systems of units.

    SYSTEM     BASIC UNIT              SYMBOL            EXAMPLE
    ============================================================

       SI      Centigrade               deg_C          32 deg_C;
       US      Fahreheit                deg_F          32 deg_F;

Example 8 : The statements

    AverageTemp = 10 deg_C;

    print "Average Daily Temperature = ", AverageTemp, "\n";
    print "Average Daily Temperature = ", AverageTemp (deg_F), "\n";

generate the output

    Average Daily Temperature =         10 deg_C 
    Average Daily Temperature =         50 deg_F

FORCES AND PRESSURE

Example 9 : The statements

    uniform_load  = 40 N/m;
    pressure_load = 20 kPa;

    print "SI Units : Uniform loading = ", uniform_load, "\n";
    print "           Presure loading = ", pressure_load, "\n";
    print "US Units : Uniform loading = ", uniform_load (lbf/ft), "\n";
    print "           Presure loading = ", pressure_load   (psi), "\n"

generate the output

    SI Units : Uniform loading =         40 N/m 
               Presure loading =         20 kPa 
    US Units : Uniform loading =      2.741 lbf/ft 
               Presure loading =      2.901 psi

ENERGY AND WORK

Note : Text in this section applies to ALADDIN 2.0 (scheduled for release in July 1997).

In engineering terms, "work" is the dot product of a force moving through a distance.

    SYSTEM     BASIC UNIT              SYMBOL            EXAMPLE
    ============================================================

       SI           Joule                 Jou           100 Jou;

Example 10 : In this example we compute the potential energy and kinetic energy of a mass-spring system. The script of input:

    distance  = 2 cm;     stiffness = 20 N/cm;
    velocity  = 1 m/sec;  mass      = 3 kg;

    print "Potential Energy = ", 1/2*stiffness*distance^2 (Jou), "\n";
    print "Kinetic Energy   = ", 1/2*mass*velocity^2      (Jou), "\n";

generate the output

    Potential Energy =        0.4 Jou 
    Kinetic Energy   =        1.5 Jou

Units of N.m are compatible with Jou.


PLANE ANGLE

Note : Text in this section applies to ALADDIN 2.0 (scheduled for release in July 1997).

Planar angles are measured in units of radians and degrees.

Example 11 : The statements

    w = 0.5*PI rad/sec;
    print "Circular Freq = ", w, "\n";
    print "Circular Freq = ", w (deg/sec), "\n";
    print "Period        = ", (2*PI rad)/w, "\n"

generate the output

    Circular Freq =      1.571 rad/sec 
    Circular Freq =         90 deg/sec
    Period        =          4 sec 

MAKING A QUANTITY DIMENSIONLESS

In the development of algorithms to solve engineering problems, sometimes it is necessary to remove the units from a physical quantity. That is why ALADDIN has a built in function QDimenLess() which removes units from quantities, as demonstrated by the following:

Example 12 : The script of code:

    print "MAKE A QUANTITY DIMENSIONLESS \n";
    x = 1 N; y = 1 cm/sec;
    z = QDimenLess(x); u = QDimenLess(y);
    print "x (with dimen)    = ", x,"\n";
    print "y (with dimen)    = ", y,"\n";
    print "x (without dimen) = ", z,"\n";
    print "y (without dimen) = ", u,"\n";

generates the output:

    MAKE A QUANTITY DIMENSIONLESS 
    x (with dimen)    =          1 N 
    y (with dimen)    =       0.01 m/sec 
    x (without dimen) =          1 
    y (without dimen) =       0.01 

ADDITION AND SUBTRACTION

Let Variable1 and Variable2 be physical quantities of the form

    Variable1 = Number1 [Length]^a . [Mass]^b . [Time]^c . [Temp]^d

    Variable2 = Number2 [Length]^e . [Mass]^f . [Time]^g . [Temp]^h

where a-g are exponent constants. The quentity sum/difference arithmetic operation

    Quantity Sum          ==>      Variable1 + Variable2 
    Quantity Difference   ==>      Variable1 - Variable2 

will be defined only if a = e, b = f, c = g , and d = h.

ALADDIN computes, by default, a physical quantity result with units matching the second physical quantity in the arithmetic operation.

Example 1 : The script of code:

    x = 1.0  m;
    y = 100 cm;

    print "x + y = ", x + y, "\n";
    print "y + x = ", y + x, "\n";

generates the output:

    x + y =        200 cm 
    y + x =          2 m 

Example 2 : An important feature of quantity operations is the check for consistent units. Suppose that we try to add a quantity with time units to a second quantity having units of length, as in the script:

    x = 1 in; y = 1 sec;
    z = x + y;

ALADDIN provides an appropriate fatal error message

    FATAL ERROR >> In Add() : Inconsistent Dimensions.
    FATAL ERROR >> Compilation Aborted.

followed by the termination of program execution.


MULTIPLICATION AND DIVISION

Let Variable1 and Variable2 be physical quantities of the form

    Variable1 = Number1 [Length]^a . [Mass]^b . [Time]^c . [Temp]^d

    Variable2 = Number2 [Length]^e . [Mass]^f . [Time]^g . [Temp]^h

where a-g are exponent constants. The quentity product and quentity division arithmetic operations are given by

    Quantity Product      ==>      Variable1 * Variable2 
    Quantity Division     ==>      Variable1 / Variable2 

The units in Variable1 * Variable2 will be

    [Length]^(a+e) . [Mass]^(b+f) . [Time]^(c+g) . [Temp]^(d+h)

and the in Variable1 / Variable2

    [Length]^(a-e) . [Mass]^(b-f) . [Time]^(c-g) . [Temp]^(d-h)

Example 3 : Suppose that the velocity and distance traveled by a free-falling object is given by

    velocity = gravity * time
    distance = 1/2 . gravity * time * time

where ``gravity'' is the acceleration due to gravity, and ``time'' is the time of free-fall.

The script of code

    gravity = 9.81 m/sec^2;
    t = 3 sec;

    print "acceleration = ", gravity , "\n";
    print "velocity     = ", gravity*t , "\n";
    print "distance     = ", 1/2*gravity*t*t, "\n";

generates the output:

    acceleration =      9.81 m/sec^2
    velocity     =      29.43 m/sec 
    distance     =      44.14 m 

Example 4 : The script of code

    force = 20    N;
    mass  =  5   kg;
    area  = 20 cm^2;

    print "force      = ", force       , "\n";
    print "mass       = ", mass        , "\n";
    print "area       = ", area (cm^2) , "\n";
    print "force/mass = ", force/mass , "\n";
    print "force/area = ", force/area , "\n";

generates the output:

    force      =         20 N 
    mass       =          5 kg 
    area       =         20 cm^2 
    force/mass =          4 m/sec^2 
    force/area =      1e+04 Pa 

EXPONENTIALS

Power operations on quantities do not work in quite the same way as the basic add, subtract, multiply and divide operations on quantities. If quantity1 and quantity2 are physical quantities, then the operation

 
      quantity1 ^ quantity2;

is defined for dimensionless quantity1 and quantity2 , and cases where one of the two operands, but not both operands, has units.

Order of Evaluation

Unlike the basic arithmetic operations, which are evaluated left-to-right, expressions involving power operations are evaluated from right-to-left. Power operations also have higher precedence than the add/subtract operators, and the muldiply/divide operators.

Example 5 : The statement

    print "2^2^3 = ", 2^2^3 , "\n"; 
generates the output
    2^2^3 =        256 
Becasue expressions involving power operations are evaluated from right-to-left, the expression
   2^2^3 
is evaluated as if it were written
   2^(2^3) 
The expression evaluates to 2^8 , and then 256 .

Example 6 : The following statements show how power operations apply to quantities having units -- the script:

    print "(2 m)^2    = ", (2 m)^2      , "\n"; 
    print "2^2 m      = ", 2^2 m        , "\n"; 
    print "10^2^2 N/m = ", 10^2^2  N/m  , "\n"; 
generates the output
    (2 m)^2    =          4 m^2 
    2^2 m      =          4 m 
    10^2^2 N/m =      1e+04 N/m 

You should should carefully note how parentheses are used to bind a unit to a constant, and how the right-to-left evaluation of numerical quantities takes precedence over the unit itself.

Example 7 : Now let's look at a family of expressions where power operations are mixed in with multiply and divide operations. The block of code

    print "0.25*2^2   m  = ", 0.25*2^2   m, "\n"; 
    print "1/2*2      m  = ",  1/2*2      m, "\n";
    print "(1/2*2     m) = ", (1/2*2     m), "\n";
    print "(1/2^2     m) = ", (1/2^2     m), "\n";
    print "(1/2^2)*(1 m) = ", (1/2^2)*(1 m), "\n";
    print "1/2^2*(1 m)   = ",  1/2^2 *(1 m), "\n";
generates the output
    0.25*2^2 m    =     1 m 
    1/2*2    m    =     1 m 
    (1/2*2  m)    =     1 m 
    (1/2^2  m)    =  0.25 1/m 
    (1/2^2)*(1 m) =  0.25 m 
    1/2^2*(1 m)   =  0.25 m 

Notice how parentheses can be used to adjust the precedence of evaluation, and how the multiply and divide operations are evaluated left-to-right.

Example 8 : In this example, we demonstrate the exponential operator on quantity variables:

    x = 100 kg; z = 10 m;
The block of input
    print "z^-3          = ", z^-3, "\n";     
    print "(x*z)^2       = ", (x*z)^2, "\n";  
generates the output
    z^-3          =  0.001 1/m^3 
    (x*z)^2       =  1e+06 m^2.kg^2 

Example 9 : Power operations such as

      x = (1 m)^(2 sec);

are illegal, and will cause ALADDIN to terminate its execution.

Physical quantities may be manipulated by some math functions ..


MATH FUNCTIONS

   FUNCTION          DESCRIPTION
   ===============================================================

   cos (x)           Compute cosine of argument x.
   sin (x)           Compute sine of argument x.
   tan (x)           Compute tangent of argument x.
   atan (x)          Compute arc-tangent of argument x.

   abs(x)            Return absolute value of quatity x.

   exp(x)            Compute exponential of quantity x.
   log(x)            Compute logarithm to "base e" of quantity x.
   log10(x)          Compute logarithm to "base 10" of quantity x.

   sqrt(x)           Return square root of quatity x.

   ===============================================================

Points to note:

Example 1 : Let's exercise a couple of the trigonometric functions. The script of code

    angle = PI/6; s = sin (angle); c = cos (angle);

    print "angle         = ", angle, "\n";
    print "sin ( angle ) = ", s, "\n";
    print "cos ( angle ) = ", c, "\n";
    print "sin^2 + cos^2 = ", s^2 + c^2, "\n";

generates the output:

    angle         =     0.5236 
    sin ( angle ) =        0.5 
    cos ( angle ) =      0.866 
    sin^2 + cos^2 =          1 
Here we compute the sine and cosine for 30 degrees, and verify that sine squared plus cos squared equals 1 (and thank goodness it does).

Example 2 : Arguments to trigonometric functions may also be defined in terms of degrees. The script of code:

    a1 = 30 deg;
    a2 = 15 deg;

    s1 = sin (a1); c1 = cos (a1);
    s2 = sin (a2); c2 = cos (a2);

    print "sin( a1 + a2 )                        = ", sin(a1 + a2), "\n";
    print "sin(a1) * cos(a2) + sin(a2) + cos(a1) = ", s1*c2 + s2*c1,"\n";

generates the output:

    sin( a1 + a2 )                        =     0.7071 
    sin(a1) * cos(a2) + sin(a2) + cos(a1) =     0.7071

and numerically verifies the well known double angular formula.

Example 3 : The abs function returns the absolute value of a physical quantity -- as a case in point, the script

    x = -10 cm/sec;
    print "x = ", x , " abs(x) = ", abs(x) , "\n";
    print "x = ", x , " abs(x) = ", abs(x) (cm/sec), "\n";

generates the output:

    x =       -0.1 m/sec  abs(x) =        0.1 m/sec 
    x =       -0.1 m/sec  abs(x) =         10 cm/sec 

You should notice how in the second example, we have cast the output of abs to an appropriate set of scaled units before printing.

Example 4 : The log and log10 functions compute the logarithms on physical quantities, after their units have been truncated. For example, the two-part script:

x = 100;
print "x = ", x , "   log(x) = ", log(x) , "\n";
print "x = ", x , " log10(x) = ", log10(x) , "\n";

y = 1000 cm^2;
print "y = ", y , "   log(y) = ", log(y) , "\n";
print "y = ", y , " log10(y) = ", log10(y) , "\n";

generates the output:

x =        100    log(x) =      4.605 
x =        100  log10(x) =          2 
y =        0.1 m^2    log(y) =     -2.303 
y =        0.1 m^2  log10(y) =         -1

The result of a logarithm computation is a dimensionless physical quantity.

Example 5 : The sqrt function returns the square root of a physical quantity, with units adjusted accordingly. As a case in point, the script

x = 16 cm^2;
print "x = ", x , " sqrt(x) = ", sqrt(x) , "\n";
print "x = ", x (cm^2), " sqrt(x) = ", sqrt(x) (cm), "\n";

generates the output:

x =     0.0016 m^2  sqrt(x) =       0.04 m 
x =         16 cm^2  sqrt(x) =          4 cm 

Now let's try to compute the square root of a negative physical quantity. If we redefine

x = - 16 cm^2;

and re-run the sqrt demonstration script, the new block of generated output is:

x =    -0.0016 m^2  sqrt(x) = sqrt: DOMAIN error
ERROR >> argument out of domain in file 'input-temp' near line 3
FATAL ERROR >> "In Warning()"

After telling you where the cause of the run-time error seems to be, ALADDIN terminates its execution.


RANDOM NUMBERS

   FUNCTION          DESCRIPTION
   ================================================================

   random ()         Return a random number selected from a uniform
                     distribution covering [0,1].

   ===============================================================

Points to note:

Example 6 : Here we simply

    print "Random 1 : ", random(), "\n";
    print "Random 2 : ", random(), "\n";
    print "Random 3 : ", random(), "\n";
    print "Random 4 : ", random(), "\n";

call random() four times, and print the results:

    Random 1 :     0.9187 
    Random 2 :     0.3461 
    Random 3 :     0.5872 
    Random 4 :     0.6017 

Because the seed for the random number generator is initialized by the "time()" function, each sequence of random numbers will be different (we should probably change this in a future version).

Logical and relational operands are used in questions that have a true/false answer. Expressions composed of logical and relational operands play a central role in ALADDIN programming constructs for the branching and looping control of problem solving procedure.


LOGICAL OPERANDS

ALADDIN supports three logical operands on physical quantities. They are:

    OPERAND           DESCRIPTION
    ====================================================

    ||                Logical "or"
    &&                Logical "and"
     !                Logical "not"

    ====================================================

Expressions involving logical operands will evaluate to either true or false. In keeping with other programming languages, "true" is represented by a non-zero number, and "false" by zero.

Example 1 : In our first example,

    print "2 && 3   : ", 2 && 3, " (true)  \n";
    print "2 && 0   : ", 2 && 0, " (false) \n";
    print "2 || 0   : ", 2 || 0, " (true)  \n";
    print "0 || 0   : ", 0 || 0, " (false) \n";
    print "!0       : ", !0    , " (true)  \n";
    print "!2       : ", !2    , " (false) \n";
we simply evaluate logical operands for a variety of non-zero and zero number combinations. The generated output is:
    2 && 3   :          1  (true)  
    2 && 0   :          0  (false) 
    2 || 0   :          1  (true)  
    0 || 0   :          0  (false) 
    !0       :          1  (true)  
    !2       :          0  (false)

Example 2 : Logical operand also work on physical quantities having units. For example, the script

    print "2 cm && 3 cm : ", 2 cm && 3 cm, " (true)  \n";
    print "0 cm/sec && 3 cm^2  : ", 0 cm && 3 cm, " (false)  \n";
generates the output
    2 cm && 3 cm :           1  (true)  
    0 cm/sec && 3 cm^2  :    0  (false)
In each case, the logical operand in applied to the numerical component of the physical quantity alone (the units play no role in the operand's evaluation).

RELATIONAL OPERANDS

The relational operands are:

    OPERAND           DESCRIPTION
    ====================================================

    ==                Test for "identical equality"
     >                Greater than
     <                Less than
    >=                Greater than or equal to
    <=                Less than or equal to

    ====================================================

Example 3 : Expressions involving relational opeators may be applied directly to physical quantities -- for example, the script:

    print "2 cm == 3 cm : ", 2 cm == 3 cm, " (false)  \n";
    print "2 cm >  3 cm : ", 2 cm >  3 cm, " (false)  \n";
    print "2 cm <  3 cm : ", 2 cm <  3 cm, " (true)   \n";
    print "2 cm <= 3 cm : ", 2 cm <= 3 cm, " (true)   \n";
generates the output
    2 cm == 3 cm :          0  (false)  
    2 cm >  3 cm :          0  (false)  
    2 cm <  3 cm :          1  (true)   
    2 cm <= 3 cm :          1  (true)

Example 4 : A key feature of ALADDIN's units package is checking of compatible units before a logical or relational operand is evaluated. In the following script, we show how sets of mixed units may be used in relational expressions -- the script of code:

    print "20 cm  == 200 mm : ", 20 cm == 200 mm, " (true)  \n";
    print "12 in   >   1 ft : ", 12 in  >   1 ft, " (false)  \n";
    print "36 in  => 100 cm : ", 36 in  > 100 cm, " (false)  \n";
generates the output
    20 cm  == 200 mm :          1  (true)  
    12 in   >   1 ft :          0  (false)  
    36 in  => 100 cm :          0  (false)  
Everything works as expected. Now suppose that we make a small error, as in the script:
    x = 20 cm;
    y = 20 cm/sec;

    print "20 cm => 20 cm/sec : ", x > y , " (huh ???)  \n";
ALADDIN generates the output
    20 cm => 20 cm/sec :
    FATAL ERROR >> "In Quantity_Gt() : Inconsistent Dimensions"
and terminates its execution.

LOGICAL AND RELATIONAL EXPRESSIONS

Of course logical and relational operands may be combined. Here are a couple of examples:

Example 5 : The script

    x = 20 cm > 10 cm;
    y = 20 sec >= 100 sec ;

    print "x : ", x , " (true)  \n";
    print "y : ", y , " (false)  \n";
    print "x || y : ", x || y , " (true)  \n";
    print "x && y : ", x && y , " (false)  \n";
generates the output
    x :          1  (true)  
    y :          0  (false)  
    x || y :          1  (true)  
    x && y :          0  (false)
Note : Further examples may be found in the sections on branching constructs and looping constructs.

The if and if-then-else statements provide branching control in ALADDIN's problem solving procedures.


THE if STATEMENT

Syntax : The if statement syntax is:

    if ( logical expression ) {
         statement 1 ;
         statement 2 ;
         ......
         statement N ;
    }

Statements 1 through N will be executed only if the logical expression evaluates to "true." Readers should also note that unlike the C programming language, the braces are required even if N equals 1.

Example 1 : The script

    x = 1 ksi;                  
    if ( x < 10 ksi ) {
         print " x = ", x ,"\n";
    }

generates the output

    x =      1 ksi

Example 2 : Logical expressions can be a composition of logical and relational operands. For example, the script:

    time = 20 sec;
  
    if( time > 0 sec && time <= 100 sec ) {
        print " time = ", time, "\n";
    }
generates the output
    time =       20 sec 

Example 3 : ALADDIN checks compatibility of units before attempting to evaluate the logical expression. The run-time execution of the script

    x = 1 ksi;                  
    if ( x < 10 cm ) {
         print " x = ", x ,"\n";
    }
fails because the units of x and 10 cm are incompatible.

THE if-then-else STATEMENT

Syntax : The if-then-else statement syntax is:

    if ( logical expression ) then {
         statement 1 ;
         ......
         statement K ;
    } else {
         statement K+1 ;
         ......
         statement N ;
    }

Statements 1 through K will be executed if the logical expression evaluates to "true." Otherwise, statements K+1 through N will be executed.

Example 4 : The script

    x = 10 ksi; y = 1 MPa;
    if ( x < 10 ksi ) then {
         print " x = ", x ,"\n";
    } else {
         print " y = ", y ,"\n";                    
    }

generates the output

    y =        1 MPa

The while and for statements provide looping control in ALADDIN's problem solving procedures. The break statement forces an early exit from a looping construct.


THE while LOOPING CONSTRUCT

Syntax : The while statement has the syntax :

    while ( boolean expression ) {
         statement 1 ;
         statement 2 ;
         ......
         statement N ;
    }

The statements 1 through N will be executed repeatedly while the boolean expression evaluates to "true."

Example 1 : The script

    while ( 1 ) {
       print " x = ", x ,"\n";
    }

contains an infinite loop that generates the output

    x =      1 ksi
    x =      1 ksi
    x =      1 ksi

    .... and on forever ....

Example 2 : Looping constructs are essential for the efficient generation of finite element meshes. Suppose, for example, that we need to generate a line of nodal coordinates. The script

    x = 0 ft;
    while( x <= 10 ft ) {
        print " x coord = ", x , "\n";
        x = x + 2 ft;
    }

generates the output

    x coord =          0 ft 
    x coord =          2 ft 
    x coord =          4 ft 
    x coord =          6 ft 
    x coord =          8 ft 
    x coord =         10 ft 

Example 3 : The coordinates in a two-dimensional finite element mesh can be generated with nested while loops. For example, the script:

    node = 0; y = 0 ft;
       while( y <= 9 ft ) {
       x = 0 ft;
       while( x <= 6 ft ) {
           print " coord ", node, "[x,y] = [", x, y, "]\n";
           node = node + 1;
           x = x + 2 ft;
       }
       y = y + 3 ft;
    }

generates 15 lines of output. An abbreviated summary is:

 coord          0 [x,y] = [         0 ft          0 ft ]
 coord          1 [x,y] = [         2 ft          0 ft ]
 coord          2 [x,y] = [         4 ft          0 ft ]
 coord          3 [x,y] = [         6 ft          0 ft ]
 coord          4 [x,y] = [         0 ft          3 ft ]

 ..... lines of output removed .....

 coord         14 [x,y] = [         4 ft          9 ft ]
 coord         15 [x,y] = [         6 ft          9 ft ]

THE for LOOPING CONSTRUCT

Syntax : The for statement has the syntax :

    for ( initial statements; boolean expression; increment statements ) {
       statement 1 ;
       statement 2 ;
       ......
       statement N ;
    }

The four steps of for-loop operation are:

Example 4 : The script of code

    for( x = 0 ft; x <= 9 ft; x = x + 3 ft ) {
         print "x = ", x, "\n";
    }
generates the output
    x =          0 ft 
    x =          3 ft 
    x =          6 ft 
    x =          9 ft 
and demonstrates behavior in a basic for loop construct.

Example 5 : The initial and increment statement blocks may contain multiple statements, separated by commas. In the script

    for( x = 0 ft, y = 0 ft; x <= 9 ft && y <= 4 ft;
         x = x + 3 ft, y = y + 2ft) {
         print "x = ", x, ": y = ", y, "\n";
    }
the initialization statements are
    x = 0 ft, y = 0 ft;
and the increment statement are
    x = x + 3 ft, y = y + 2ft;
The generated output is
   x =          0 ft : y =          0 ft 
   x =          3 ft : y =          2 ft 
   x =          6 ft : y =          4 ft
The looping construct terminates when y is incremented to 6 ft and the relational operation to y <= 6 ft evaluates to false.

Example 6 : The initialization and increment statement blocks may be empty. For example, the script of code

    x = 0 ft; y = 0 ft; 
    for( ; x <= 9 ft && y <= 4 ft; ) {
         print "x = ", x, ": y = ", y, "\n";
         x = x + 3 ft;
         y = y + 2 ft;
    }
generates the same block of output as in Example 5.

THE break STATEMENT

The break statement forces an early exit from one layer of looping constructs. We will demonstrates its behavior via a series of examples.

Example 6 : The script of statements:

    for( x = 0 ft ; x <= 9 ft ; x = x + 3 ft ) {
         print "x = ", x, "\n";
         if (x == 6 ft) 
             break;
      
    }
    print "breaking from loop \n";
    print "x = ", x, "\n";
generates the same block of output
    x =          0 ft 
    x =          3 ft 
    x =          6 ft 
    breaking from loop 
    x =          6 ft
and shows how the program execution exits the looping construct before the boolean expression has an opportunity to evaluate to false

Example 7 : What about 2 levels of looping construct ? In the script of statements:

    for( x = 0 ft ; x <= 6 ft ; x = x + 3 ft ) {
        for( y = 0 ft ; y <= 9 ft ; y = y + 3 ft ) {
             print "x = ", x, "y = ", y, "\n";
             if (y == 6 ft) 
                 break;
        }
    }
    print "finished looping construct\n";
    print "x = ", x, "\n";
the output is
    x =          0 ft y =          0 ft 
    x =          0 ft y =          3 ft 
    x =          0 ft y =          6 ft 
    x =          3 ft y =          0 ft 
    x =          3 ft y =          3 ft 
    x =          3 ft y =          6 ft 
    x =          6 ft y =          0 ft 
    x =          6 ft y =          3 ft 
    x =          6 ft y =          6 ft 
    finished looping construct
    x =          9 ft
You should notice how the values of y in the inner for loop never reach 9 ft because the relational operation
    y == 6 ft 
always evaluates to "true" beforehand. The break statement forces the program execution to exit one level of looing constructs.

MATRICES OF NUMBERS

In ALADDIN, floating point numbers should be thought of as dimensionless physical quantities.

Example 1 : Small matrices of numbers may be defined as follows :

   /* [a] : Define a (1x3) matrix called X */

   X = [ 1, 2, 3 ];

   /* [b] : Define a (3x1) matrix called Y */

   Y = [ 1; 2; 3 ];

   /* [c] : Define a (2x2) matrix called Z */

   Z = [ 1, 3; 2, 4];

Each matrix definition begins with a "[" and ends with a "]". The elements along a matrix row are separated by commas. A semi-colon between matrix elements indicates the end of a matrix row. Don't forget to terminate the matrix definition with a semi-colon !

The name of the matrix is given on the left-hand side of the assignment (i.e. "=") operator. Please ensure that you choose a matrix name that will not conflict with a keyword or function name in the ALADDIN language.

Matrices may be defined over multiple lines of input.

Example 2 : Matrix Z could have also been written:

   /* [c] : Define a (2x2) matrix Z */

   Z = [ 1, 3;
         2, 4 ];

MATRIX ELEMENTS

To extract an element of a matrix, simply append the matrix name with the row and column number enclosed in brackets.

Example 3 : From example 2, we can write

    print "Z[1][1] = ", Z[1][1], "\n";
    print "Z[2][1] = ", Z[2][1], "\n";
    print "Z[1][2] = ", Z[1][2], "\n";
    print "Z[2][2] = ", Z[2][2], "\n";
and the output will be
    Z[1][1] =       1
    Z[2][1] =       2
    Z[1][2] =       3
    Z[2][2] =       4

In ALADDIN, matrix row and column numbers begin with one. The matrix elements,

    Z[1][1], Z[1][2] ... 

are physical quantities. If a matrix element has units (see examples below) then the units will be printed along with the numerical value of the matrix element.


MATRICES OF QUANTITIES WITH UNITS

ALADDIN supports the specification of rectangular arrays of physical quantities.

Example 4 : The statement:

    box = [ 3 m, 6 m, 9 m ]; 

defines a (1x3) matrix containing length quantities. Now let's suppose that the elements of Z represent the side lengths of a rectangular box. The box volume can be computed and printed by simply typing:

    print "volume of box = ", box[1][1]*box[1][2]*box[1][3], "\n";

The output is:

    volume of box =       162 m^3

PRINTING MATRICES

The PrintMatrix() function can print one or more matrices.

Example 5 : The pair of statements:

    PrintMatrix(X);
    PrintMatrix(Y,Z);

print matrices X, Y and Z, as defined in Examples 1 and 2. The output is:

    MATRIX : "X"

    row/col                  1            2            3   
            units                                          
       1            1.00000e+00  2.00000e+00  3.00000e+00

    MATRIX : "Y"

    row/col                  1   
            units                
       1            1.00000e+00
       2            2.00000e+00
       3            3.00000e+00

    MATRIX : "Z"

    row/col                  1            2   
            units                             
       1            1.00000e+00  3.00000e+00
       2            2.00000e+00  4.00000e+00

In the second function call, matrices Y and Z are printed with only one function call to PrintMatrix().

Example 6 : The function call (from Example 3)

    PrintMatrix( box );

generates the output:

    MATRIX : "box"

    row/col                  1            2            3   
            units            m            m            m   
       1            3.00000e+00  6.00000e+00  9.00000e+00

Notice that the units of "m" are aligned along each column.


FUNCTIONS FOR MATRIX ALLOCATION

The explicit definition of matrices is really only practical for matrix sizes less than about 3-5 rows and columns. For all other cases, ALADDIN provides a suite of functions for the generation of matrices commonly used in numerical, matrix, and/or finite element applications.

Functions for the Defintion of Matrices

    FUNCTION		PURPOSE
    =========================================================================

    Matrix([s,t])        Allocate a s x t matrix
    Diag([s, n])         Allocate s x s diagonal matrix with n along diagonal			
    One([s, t])          Allocate s x t matrix full of ones
    One([s])             Allocate s x s matrix full of ones
    Zero([s])            Allocate s x s matrix full of zeros
    Zero([s, t])         Allocate s x t matrix full of zeros

Example 7 : The following script exercises each of these matrix allocation functions:

    /* [a] : Allocate a 20 by 30 matrix */

    W = Matrix([20,30]);

    /* [b] : Allocate a 1 by 30 matrix full of zeros */

    X = Zero([1,30]);

    /* [c] : Allocate a 30 by 30 matrix full of zeros */

    X = Zero([30,30]);
    Y = Zero([30]);

    /* [d] : Allocate a matrix full of zeros, the size the same as [W] */

    X = Zero( Dimension(W) );

    /* [e] : Allocate a 30 by 30 matrix full of ones */

    Y = One([30]);
    X = One([30,30]);

    /* [f] : Allocate a 30 by 30 diagonal matrix with 2 along the */
    /*	    the diagonal and a 44 by 44 identity matrix          */

    X = Diag([30,2]);
    Y = Diag([44,1]); 

ROW UNITS AND COLUMN UNITS

The functions

    FUNCTION		PURPOSE
    ======================================================

    ColumnUnits(A, [u])  Assign column units u to matrix A
    RowUnits(A, [u])     Assign row units u to matrix 

allow for the definition of matrices with specified row and column units.

Example 8 : The block of statements:

    stiff = Matrix([3,3]);
    stiff = ColumnUnits ( stiff, [ N/m ] );

    stiff [1][1] = 10 N/m;
    stiff [2][2] = 20 N/m;
    stiff [3][3] = 30 N/m;

    PrintMatrix(stiff);

generates the output:

    MATRIX : "stiff"

    row/col                  1            2            3   
            units          N/m          N/m          N/m   
       1            1.00000e+01  0.00000e+00  0.00000e+00
       2            0.00000e+00  2.00000e+01  0.00000e+00
       3            0.00000e+00  0.00000e+00  3.00000e+01

The Matrix() function takes care of the dynamic allocation of memory for "stiff" ... and the ColumnUnits() function assigns the units of Newtons per meter to each column.

Example 9 : Only one command line is needed to specify a matrix having a common set of column units. As a case in point, our test matrix "stiff" could have been allocated with :

    stiff = ColumnUnits ( Matrix([3,3]) , [ N/m ] );

PRINTING MATRICES WITH SCALED ROW/COLUMN UNITS

Note : Text in this section applies to ALADDIN 2.0 (scheduled for release in July-August 1997).

Matrices may be printed with either the row units (or column units) scaled to a desirable range.

The function syntax is:

    PrintMatrixCast ( matrix1, ..., matrixN, [ units matrix ] );

Here matrix1 through matrixN are N matrices having the row or column units listed in

    [ units matrix ]

Column units are scaled by supplying a units matrix having one row and more than one column. Row units are scaled by supplying a units matrix having more than one row and only one column.

An error condition will occur if terms in the units matrix do not match any of the terms in matrices matrix1 through matrixN.

Example 10 : The block of statements:

    stiff = Matrix([3,3]);
    stiff = ColumnUnits ( stiff, [ N/m , N/m , N/rad ] );

    stiff [1][1] = 10 N/m;
    stiff [2][2] = 20 N/m;
    stiff [3][3] = 30 N/rad;

    PrintMatrixCast(stiff, [ N/cm, kN/rad ] );

generates the output:

    MATRIX : "stiff"

    row/col                  1            2            3   
             units         N/cm         N/cm       kN/rad   
       1            1.00000e-01  0.00000e+00  0.00000e+00
       2            0.00000e+00  2.00000e-01  0.00000e+00
       3            0.00000e+00  0.00000e+00  3.00000e-02

Example 11 : The block of statements:

    spectra = Matrix( [3,2] );
    spectra = ColumnUnits( spectra, [sec],      [1]);
    spectra = ColumnUnits( spectra, [cm/sec^2], [2]);

    spectra [ 1][1] = 0.0 sec;  spectra [ 1][2] = 981.0*0.15 cm/sec/sec;
    spectra [ 2][1] = 0.1 sec;  spectra [ 2][2] = 981.0*0.18 cm/sec/sec;
    spectra [ 3][1] = 0.2 sec;  spectra [ 3][2] = 981.0*0.25 cm/sec/sec;

    PrintMatrix (spectra);

    spectra = ColumnUnits( spectra, [ft/sec^2], [2]);
    PrintMatrixCast (spectra , [sec, cm/sec^2] );
    PrintMatrixCast (spectra , [sec, in/sec^2] );

generates the output:

    MATRIX : "spectra"

    row/col                  1            2   
            units          sec      m/sec^2   
       1            0.00000e+00  1.47150e+00
       2            1.00000e-01  1.76580e+00
       3            2.00000e-01  2.45250e+00

    MATRIX : "spectra"

    row/col                  1            2   
            units          sec     cm/sec^2   
       1            0.00000e+00  1.47150e+02
       2            1.00000e-01  1.76580e+02
       3            2.00000e-01  2.45250e+02

    MATRIX : "spectra"

    row/col                  1            2   
            units          sec     in/sec^2   
       1            0.00000e+00  5.79331e+01
       2            1.00000e-01  6.95197e+01
       3            2.00000e-01  9.65551e+01

The first and second columns of "spectra" contain quantities of time and acceleration, respectively. Here we define the accelerations in terms of cm/sec/sec, and then use the PrintMatrixCast function to scale the output to an appropriate set of units.

For more details and examples, please consult Technical Report T.R. 95-74.


ADDITION AND SUBTRACTION

Let X be a (mxn) matrix and Y be a (pxq) matrix. The matrix sum "X + Y" is defined when:

Example 1 : In the following script we define two (2x2) matrices, compute their sum, and print all of the participating matrices. The script:

    /* [a] : Define (2x2) matrices "A" and "B" */ 

    A = [ 1, 2; 3, 4];
    B = [ 3, 4; 5, 6];

    /* [b] : Compute and print the matrix sum "A+B" */ 

    C = A+B;
    PrintMatrix( A, B, C);

generates the output:

    MATRIX : "A"

    row/col                  1            2
            units
       1            1.00000e+00  2.00000e+00
       2            3.00000e+00  4.00000e+00

    MATRIX : "B"

    row/col                  1            2
            units
       1            3.00000e+00  4.00000e+00
       2            5.00000e+00  6.00000e+00

    MATRIX : "C"

    row/col                  1            2
            units
       1            4.00000e+00  6.00000e+00
       2            8.00000e+00  1.00000e+01

Example 2 : The script of code:

    /* [a] : Initialize (4x1) matrix "X" */

    X = Diag([4, 1]);                    
    X = ColumnUnits(X, [ksi, ft,  N,  m]);
    X =    RowUnits(X, [psi, in, kN, mm]);

    /* [b] : Initialize (4x1) matrix "Y" */

    Y = One([4]);   
    Y = ColumnUnits(Y, [psi, in, kN, km]);
    Y =    RowUnits(Y, [ksi, ft,  N, mm]);

    /* [c] : Compute and print matrix sum "Z = X + Y" */

    Z = X + Y;
    PrintMatrix(X, Y, Z);

generates the output:

    MATRIX : "X"

    row/col                  1            2            3            4   
            units          ksi           ft            N            m   
       1      psi   1.00000e+00  0.00000e+00  0.00000e+00  0.00000e+00
       2       in   0.00000e+00  1.00000e+00  0.00000e+00  0.00000e+00
       3       kN   0.00000e+00  0.00000e+00  1.00000e+00  0.00000e+00
       4       mm   0.00000e+00  0.00000e+00  0.00000e+00  1.00000e+00

    MATRIX : "Y"

    row/col                  1            2            3            4   
            units          psi           in           kN           km   
       1      ksi   1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
       2       ft   1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
       3        N   1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
       4       mm   1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00

    MATRIX : "Z"

    row/col                  1            2            3            4   
            units          ksi           ft            N            m   
       1      psi   2.00000e+00  8.33333e+01  1.00000e+06  1.00000e+06
       2       in   1.20000e-02  2.00000e+00  1.20000e+04  1.20000e+04
       3       kN   1.00000e-06  8.33333e-05  2.00000e+00  1.00000e+00
       4       mm   1.00000e-03  8.33333e-02  1.00000e+03  1.00100e+03

Matrix addition is commutative (i.e. X + Y = Y + X).


MULTIPLICATION

Let X be a (mxn) matrix and Y be a (pxq) matrix. The matrix product X.Y is defined when:

Example 3 : Here is an example of multiplication of matrices containing dimensionless physical quantities -- the script of input:

    /* [a] : Define (3x3) X matrix, and (3x4) Z matrix */

    X = [ 2, 3, 5;
          4, 6, 7;
         10, 2, 3];
    Z = One([3, 4]);

    /* [b] : Compute U = X*Z and print all matrices */

    U = X*Z;
    PrintMatrix(X, Z, U);

The generated output is:

    MATRIX : "X"

    row/col                  1            2            3   
            units                                          
       1            2.00000e+00  3.00000e+00  5.00000e+00
       2            4.00000e+00  6.00000e+00  7.00000e+00
       3            1.00000e+01  2.00000e+00  3.00000e+00

    MATRIX : "Z"

    row/col                  1            2            3            4   
        units                                                       
       1            1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
       2            1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00
       3            1.00000e+00  1.00000e+00  1.00000e+00  1.00000e+00

    MATRIX : "U"

    row/col                  1            2            3            4   
            units                                                       
       1            1.00000e+01  1.00000e+01  1.00000e+01  1.00000e+01
       2            1.70000e+01  1.70000e+01  1.70000e+01  1.70000e+01
       3            1.50000e+01  1.50000e+01  1.50000e+01  1.50000e+01

Example 4 : The following script demonstrates multiplication of matrices containing physical quantities that have dimensions:

    /* [a] : Define (3x3) stiffness matrix */

    stiff = [ 2, 3, 5;     
              4, 6, 7;
             10, 2, 3];
    stiff = ColumnUnits( stiff, [N/m, N/m, N/rad ]);
    stiff =    RowUnits( stiff, [m] , [3]);

    /* [b] : Define (3x1) displacement matrix */

    displ = One([3, 1]);   
    displ = RowUnits( displ, [m, m, rad]);

    /* [c] : Compute (3x1) external load vector/matrix */

    load = stiff*displ; 

    PrintMatrix( stiff, displ, load );

The generated output is:

    MATRIX : "stiff"

    row/col                  1            2            3   
            units          N/m          N/m        N/rad   
       1            2.00000e+00  3.00000e+00  5.00000e+00
       2            4.00000e+00  6.00000e+00  7.00000e+00
       3        m   1.00000e+01  2.00000e+00  3.00000e+00

    MATRIX : "displ"

    row/col                  1   
            units                
       1        m   1.00000e+00
       2        m   1.00000e+00
       3      rad   1.00000e+00

    MATRIX : "load"

    row/col                  1   
            units                
       1        N   1.00000e+01
       2        N   1.70000e+01
       3      N.m   1.50000e+01

SCALING A MATRIX BY A QUANTITY

Let A be a matrix and q be a physical quantity. ALADDIN's supports pre- and post-multiplication of A by quantity q (i.e. q.A and A.q), and division of A by a quantity (i.e. A/q).

Example 5 : The script

    /* [a] : Define (3x2) matrix and velocity quantity */

    coord = [ 1, 2;
              3, 4;
              5, 6 ];
    velocity = 2 in/sec;

    /* [b] : Compute and print the matrix sum */

    Gv = coord*velocity;
    PrintMatrix( coord, Gv );

generates the output

    MATRIX : "coord"

    row/col                  1            2   
            units                             
       1            1.00000e+00  2.00000e+00
       2            3.00000e+00  4.00000e+00
       3            5.00000e+00  6.00000e+00

    MATRIX : "Gv"

    row/col                  1            2   
            units       in/sec       in/sec   
       1            2.00000e+00  4.00000e+00
       2            6.00000e+00  8.00000e+00
       3            1.00000e+01  1.20000e+01

and demonstrates post-multiplication of a matrix by a quantity. Notice how the units of "velocity" are transformed into the matrix-quantity product.

Similar results are obtained for matrix-quantity expressions

       velocity*coord;  coord/velocity;
and so forth.

For more details and examples, please consult Technical Report T.R. 95-74.


DIMENSIONS OF A MATRIX

The number of rows and columns in a matrix may be extracted, and used in the construction of general purpose problem solving procedures. For example, a common use for matrix dimensions is in setting the beginning and end points for looping constructs in numerical algorithms.

    FUNCTION           PURPOSE
    =====================================================================

    Dimension( A )     Extract the number of rows and columns in matrix
                       "A" and return the result in a (1x2) matrix.

                       Element [1][1] holds the number of rows in "A"
                       Element [1][2] holds the number of columns in "A"

Example 1 : The script of code:

    Z = One([ 13, 20]);    /* Allocate 3 by 20 matrix called Z   */
    size = Dimension(Z);   /* Extract and print dimensions of Z  */
   
    print "Rows in [Z]     = ",size[1][1] ,"\n"; 
    print "Columns in [Z]  = ",size[1][2] ,"\n";

generates the output:

    Rows in [Z]    = 1.3000e+01
    Columns in [Z] = 2.0000e+01

MATRIX COPY

Let X be a (mxn) matrix. A copy of X can be made by simply writing

    Y = X;

There are situations, however, where a matrix copy is needed, but without an assignment. The matrix function

    FUNCTION                  DESCRIPTION
    ================================================================
    Copy(A)                   Takes one matrix argument, and returns 
                              a matrix copy.

Of course, Copy() can be used with an assignment.

Example 2 : The script of code:

    /* [a] : Define and print (1x3) test matrix */ 

    response = [ 1 sec, 2 cm/sec, 3 cm/sec^2 ];
    PrintMatrix( response );

    /* [b] : Copy and print "response" matrices */ 

    copy1 = response;
    copy2 = Copy ( response );
    PrintMatrix( copy1, copy2 );

generates the output:

    MATRIX : "response"

    row/col                  1            2            3   
            units          sec        m/sec      m/sec^2   
       1            1.00000e+00  2.00000e-02  3.00000e-02

    MATRIX : "copy1"

    row/col                  1            2            3   
            units          sec        m/sec      m/sec^2   
       1            1.00000e+00  2.00000e-02  3.00000e-02

    MATRIX : "copy2"
   
    row/col                  1            2            3   
            units          sec        m/sec      m/sec^2   
       1            1.00000e+00  2.00000e-02  3.00000e-02

Example 3 : The script of code:

    /* [a] : Define and print (1x3) test matrix */ 

    response = [ 1 sec, 2 cm/sec, 3 cm/sec^2 ];
    PrintMatrix( response );

    /* [b] : Copy and print "response" matrices */ 

    copy1 = response;
    PrintMatrix( copy1, Copy (response) );

generates the same block of output as in Example 2.


MATRIX TRANSPOSE

Let "A" be a (mxn) matrix. The matrix transpose of "A" is a (nxm) matrix with the rows and columns of A interchanged.

    FUNCTION          DESCRIPTION
    ================================================================
    Trans(A)          Generates a new matrix corresponding to the
                      matrix transpose of "A"

Example 4 : The script of code:

    /* [a] : Define and print (1x3) test matrix */ 

    response = [ 0 sec, 0 cm/sec, 0 cm/sec^2 ;
                 1 sec, 2 cm/sec, 3 cm/sec^2 ];
    PrintMatrix( response );

    /* [b] : Compute and print transpose of "response" */ 

    transpose1 = Trans(response);
    PrintMatrix( transpose1 );

generates the output:

    MATRIX : "response"

    row/col                  1            2            3   
            units          sec        m/sec      m/sec^2   
       1            0.00000e+00  0.00000e+00  0.00000e+00
       2            1.00000e+00  2.00000e-02  3.00000e-02

    MATRIX : "transpose1"

    row/col                  1            2   
            units                             
       1      sec   0.00000e+00  1.00000e+00
       2    m/sec   0.00000e+00  2.00000e-02
       3  m/sec^2   0.00000e+00  3.00000e-02

MAXIMUM AND MINIMUM MATRIX ELEMENTS

    FUNCTION           PURPOSE
    ===============================================================

    Min ( A )          Return a (1x1) matrix containing the minimum 
                       matrix element in matrix "A".

    Max ( A )          Return a (1x1) matrix containing the maximum
                       matrix element in matrix "A".

Example 5 : The script

    A = [ 3.78,  9.7, -4.7,  10.50  ;
          0.00, -5.8,  0.2,  -9.34] ;

    MaxValue = Max( A );
    MinValue = Min( A );

    PrintMatrix(A);

    print "\n";
    print "Max(A) =", MaxValue ,"\n";
    print "Min(A) =", MinValue ,"\n";
generates the output
    MATRIX : "A"

    row/col                  1            2            3            4
           units
       1            3.78000e+00  9.70000e+00 -4.70000e+00  1.05000e+01
       2            0.00000e+00 -5.80000e+00  2.00000e-01 -9.34000e+00

    Max(A) =      10.5
    Min(A) =     -9.34

Note. In Aladdin 2, Max() and Min() truncate the units from the matrix element. In other words, the max/min test applies to the max/min numerical value of the matrix element. We need to change the program so that when all of the elements have the same units type, Min()/Max() also returns the units.


EUCLIDEAN NORM

The L2 norm of a row matrix or column matrix is simply the square root of the sum of the matrix elements squared.

    FUNCTION             PURPOSE
    ===================================================================

    L2Norm( A )          Compute L2 norm of either a (1xn) matrix or a
                         (nx1) matrix.

Example 6 : The script

    /* [a] : Define (1x4) test vector and compute L2 norm */ 

    testVector = [ 1, 2, 3, 4 ];
    norm  = L2Norm( testVector );

    /* [b] : Print "testVector" and L2 norm */ 

    PrintMatrix( testVector );
    print "\n";
    print "L2 norm of testVector is :", norm, "\n";
generates the output:
    MATRIX : "testVector"

    row/col                  1            2            3            4   
            units                                                       
       1            1.00000e+00  2.00000e+00  3.00000e+00  4.00000e+00

    L2 norm of testVector is :     5.477

In this page we demonstrate use of the matrix functions for the solution of linear systems of equations. We will assume readers are familiar with:

A review of the theoretical concepts needed to solve linear matrix equations may be found in Technical Report T.R. 95-74.


SOLVING LINEAR EQUATIONS [A].[x] = [b]

Let "A" be a (nxn) matrix, and "b" be a (nx1) matrix. The most straightforward way of solving the matrix equations

    [ A ].[ x ] = [ b ]
is with
    FUNCTION           PURPOSE
    ==============================================================

    Solve( A, b )      Solve matrix equations [A].[x] = [b].

    ==============================================================

Example 1 : The script of code

    /* [a] : Setup matrices [A] and [b] */ 

    A = [ 1 , 2; 3 , 4 ];
    B = [ 5 ; 11 ];

    PrintMatrix( A, B );

    /* [b] : Solve matrices and print solution vector */ 

    X = Solve (A, B);
    PrintMatrix( X );
generates the output
    MATRIX : "A"

    row/col                  1            2   
            units                             
       1            1.00000e+00  2.00000e+00
       2            3.00000e+00  4.00000e+00

    MATRIX : "B"

    row/col                  1   
            units                
       1            5.00000e+00
       2            1.10000e+01

    MATRIX : "X"

    row/col                  1   
           units                
       1            1.00000e+00
       2            2.00000e+00

Example 2 : Now let's solve a system of equations having units. The script of code (this example only works for ALADDIN 2.0, scheduled for release in July 1997):

    /* [a] : Structure Calculation */

    E = 200 MPa;
    I = 0.4 m^4;
    A = 0.1 m^2;
    L = 5.0 m;

    /* [b] : Setup and print 3x3 stiffness matrix */

    stiff = Matrix([3,3]);
    stiff = ColumnUnits( stiff, [N/m, N/m, N/rad]);
    stiff = RowUnits( stiff,     [m], [3] );

    stiff[1][1] =  E*A/L;
    stiff[2][2] = 12*E*I/L^3;
    stiff[2][3] = -6*E*I/L/L;
    stiff[3][2] = -6*E*I/L/L;
    stiff[3][3] =  4*E*I/L;

    PrintMatrix(stiff);

    /* [c] : Setup and print 3x1 force vector  */

    force = [ 0.0 N ; -5 N; 0.0 N*m ];
    PrintMatrix(force);

    /* [d] : Compute and print displacement matrix */

    displ = Solve(stiff, force ); 
    PrintMatrix(displ);
    PrintMatrixCast(displ, [ cm; deg ] );

generates the output:

    MATRIX : "stiff"

    row/col                  1            2            3   
            units          N/m          N/m        N/rad   
       1            4.00000e+06  0.00000e+00  0.00000e+00
       2            0.00000e+00  7.68000e+06 -1.92000e+07
       3        m   0.00000e+00 -1.92000e+07  6.40000e+07

    MATRIX : "force"

    row/col                  1   
            units                
       1        N   0.00000e+00
       2        N  -5.00000e+00
       3      N.m   0.00000e+00
    
    MATRIX : "displ"

    row/col                  1   
            units                
       1        m   0.00000e+00
       2        m  -2.60417e-06
       3      rad  -7.81250e-07

    MATRIX : "displ"

    row/col                  1   
        units                
       1       cm   0.00000e+00
       2       cm  -2.60417e-04
       3      deg  -4.47623e-05

In the final block of matrix output, radians are scaled to degrees, and translational displacements are expressed in "cm."


[L][U] DECOMPOSITION

In the method of LU decomposition, matrix "A" is decomposed into a product of lower and upper triangular matrices.

    Step 1 : Decompose [ A ]  ====>  [ L ].[ U ]

The solution vector "x" is computed via forward substitution and then backward susbstitution

    Step 2 : Solve [ L ].[ z ] = [ b ] : Forward  Substitution
    Step 3 : Solve [ U ].[ x ] = [ x ] : Backward Substitution

Matrix functions are provided for the LU decomposition and forward/backward substitution.

    FUNCTION             PURPOSE
    ==============================================================

    Decompose( A )       Decompose "A" into a product of lower and
                         upper triangular matrices.

    Substitution ( LU )  Compute forward/backward substitution on
                         decomposed LU matrix product.

    ==============================================================

Example 3 : The script of ALADDIN commands

    /* [a] : Setup matrices [A] and [b] */ 

    A = [ 1 , 2; 3 , 4 ];
    B = [ 5 ; 11 ];

    PrintMatrix( A, B );

    /* [b] : Solve matrices and print solution vector */ 

    LU = Decompose (A);
    X  = Substitution ( LU, B);
    PrintMatrix( X );

produces output identical to Example 1.


SOLVING [A].[x] = [b] for multiple [b]'s

In the method of LU decomposition, solutions to multiple right-hand vectors can be computed once the matrix [A] has been decomposed.

Example 4 : The script

    /* [a] : Setup matrices [A], [B1], and [B2] */

    A = [ 1 , 2; 3 , 4 ];
    B1 = [ 5 ; 11 ];
    B2 = [ 1 ;  2 ];

    PrintMatrix( A );

    /* [b] : Decompose "A" into the "LU" matrix product */

    LU = Decompose (A);

    /* [c] : Solve matrix equations via forward/backward substitution */

    X1  = Substitution ( LU, B1);
    PrintMatrix( B1, X1 );

    X2  = Substitution ( LU, B2);
    PrintMatrix( B2, X2 );
solves two sets of equations A.x = B, with A as defined in Examples 1 and 2. The abbreviated output is:
    SOLVING [A].[X1] = [B1]              SOLUTION VECTOR
    ================================================================

    MATRIX : "B1"                        MATRIX : "X1"

    row/col                  1           row/col                  1   
            units                                units                
       1            5.00000e+00             1            1.00000e+00
       2            1.10000e+01             2            2.00000e+00

    SOLVING [A].[X2] = [B2]              SOLUTION VECTOR
    ================================================================

    MATRIX : "B2"                        MATRIX : "X2"

    row/col                  1           row/col                  1   
            units                                units                
       1            1.00000e+00             1            0.00000e+00
       2            2.00000e+00             2            5.00000e-01

INVERSE OF MATRIX [A]

The inverse of a square matrix "A" is computed with

    FUNCTION           PURPOSE
    ==============================================================

    Inverse( A )       Compute inverse of square matrix [A]

    ==============================================================

Example 5 : The script of code

    /* [a] : Setup matrix [A] and compute [A]^-1 */ 

    A = [ 1,  2,  3;
          7, 10,  3;
          3,  7,  8 ];

    Ainv = Inverse (A);
    PrintMatrix( A , Ainv);

    /* [b] : Check results by computing [A].[Ainv] */ 

    PrintMatrix( A*Ainv );

defines a (3x3) matrix [A], and computes its inverse. We check that [A].[Ainv] equals [I], a (3x3) identity matrix. The output is:

    MATRIX : "A"

    row/col                  1            2            3   
            units                                          
       1            1.00000e+00  2.00000e+00  3.00000e+00
       2            7.00000e+00  1.00000e+01  3.00000e+00
       3            3.00000e+00  7.00000e+00  8.00000e+00

    MATRIX : "Ainv"

    row/col                  1            2            3   
            units                                          
       1            2.68182e+00  2.27273e-01 -1.09091e+00
       2           -2.13636e+00 -4.54545e-02  8.18182e-01
       3            8.63636e-01 -4.54545e-02 -1.81818e-01

    MATRIX : "UNTITLED"

    row/col                  1            2            3   
            units                                          
       1            1.00000e+00  0.00000e+00  0.00000e+00
       2            8.88178e-16  1.00000e+00  4.44089e-16
       3           -8.88178e-16  0.00000e+00  1.00000e+00

Example 6 : Now let's add units to matrix X in Example 5.

If [A] is defined as

    /* [a] : Setup matrix [A] and compute [A]^-1 */ 

    A = ColumnUnits ([ 1,  2,  3;
                       7, 10,  3;
                       3,  7,  8 ], [ N/m, N/m, m ]);

then the inverse of A is

    MATRIX : "Ainv"

    row/col                  1            2            3   
            units                                          
       1      m/N   2.68182e+00  2.27273e-01 -1.09091e+00
       2      m/N  -2.13636e+00 -4.54545e-02  8.18182e-01
       3      1/m   8.63636e-01 -4.54545e-02 -1.81818e-01

GENERALISED SYMMETRIC EIGENPROBLEM

In this page we demonstrate we employ ALADDIN's matrix functions to compute solutions to the generalised symmetric eigenvalue problem. More precisely, we seek solutions to the matrix equations:

    [ K ].[ phi ] = [ M ].[ phi ].[ lambda ]

    (nxn) . (nxp) = (nxn) . (nxp) . (pxp)     <=== matrix dimensions.
where
    [ K ]      = (nxn) positive definite symmetric matrix;
    [ M ]      = (nxn) positive definite symmetric matrix;
    [ phi ]    = (nxp) matrix containing 'p' eigenvectors;
    [ lambda ] = (pxp) diagonal matrix containing 'p' eigenvalues;

and p <= n . For structural dynamics applications, [K] and [M] correspond to the stiffness and mass matrices, respectively.

We will assume readers are familiar with:

A review of the theoretical concepts needed to solve linear matrix equations may be found in Bathe's Finite Element Text.


EIGENVALUES AND EIGENVECTORS

ALADDIN provides three functions for computing solutions to the generalised symmetric eigen problem:

    FUNCTION               PURPOSE
    ==================================================================

    Eigen( K, M, [p] )     Compute 'p' eigenvalues and eigenvectors.

                           The third argument is a (1x1) matrix cont-
                           aining the number of eigenvalues/vectors
                           to be computed (i.e. `p`).

                           The result is assigned to a ((n+1) x p)
                           matrix called [ eigen ].

    Eigenvalue( eigen )    Extract eigenvalues from array [ eigen ].

    Eigenvector( eigen )   Extract eigenvectors from array [ eigen ].

    ==================================================================

Example 1 : The script of code

    /* [a] : Setup (2x2) [K] matrix and [M] = [I] */ 

    K = [  2,  -2;
          -2,   4 ];

    M = [ 1,  0;
          0,  1 ];

    /* [b] : Compute eigenvalues and eigenvectors */ 

    eigen = Eigen ( K, M, [2]);
    eigenvalue  = Eigenvalue ( eigen );
    eigenvector = Eigenvector ( eigen );

    /* [c] : Print eigenvalues and eigenvectors */ 

    PrintMatrix( eigenvalue );
    PrintMatrix( eigenvector );

generates the output

    *** SUBSPACE ITERATION CONVERGED IN  2 ITERATIONS 

    MATRIX : "eigenvalue"

    row/col                  1   
           units                
       1            7.63932e-01
       2            5.23607e+00

   MATRIX : "eigenvector"

   row/col                  1            2   
           units                             
      1            1.00000e+00 -6.18034e-01
      2            6.18034e-01  1.00000e+00

In thie example we compute both eigenvalues and eigenvectors. Only 2 cycles of subspace iteration are needed to reach convergece.

Note : By writing down the characteristic equations for [K] it is and easy hand calculation to show that the eigenvalues are given by solutions to

    P (x) = (2 - x).(4 - x) - 4 = 0. ==> x_1 = 3 - sqrt(5) = 0.764.
                                         x_2 = 3 + sqrt(5) = 5.236.

References
  1. Aladdin Web Site: See http://www.isr.umd.edu/~austin/aladdin.html, 2003.