OVERVIEW

Mathnium is a program for interactive numerical computations. With its comprehensive library of functions for a variety of problems in applied mathematics, and with its facilities for the definition and manipulation of arrays of numbers as basic data objects, Mathnium allows you to solve numerical problems rather painlessly and without a great deal of programming effort.

You can also extend the capabilities of Mathnium by using it as a programming language. While it is not a general purpose language, its facilities are more than sufficient for you to be able to implement your algorithms in a manner which is not only extremely concise but also very close to the mathematical description of the algorithms.

Mathnium is written in Java, and it provides convenient language constructs which do not require you to actually write formal Java programs to uses existing Java classes. As a result, very large number of existing libraries of Java classes, most of them available from open source projects, can be very easily used within the Mathnium environment.

Mathnium is inspired by Matlab, a very successful and versatile language for scientific computing. Although it should be possible to execute within the Mathnium environment a Matlab program that does not use any functions that are not in the Mathnium library, no attempt has been made (or is planned) to make Mathnium completely compatibile with Matlab.

INTERACTING WITH MATHNIUM

The simplest use of Mathnium is as a calculator: you type in a expression or a group of expressions separated by commas, and Mathnium types the results back.
>>2+3, 3^2                # A simple statement 
      5     9 
>>s=0 
>>for i=1:10          # A control structure 
>    s=s+i 
>  end 
>>s 
     55  
It is also possible to execute a set of commands from a text file by using the function exec. If you issue the command

exec("cfile") 

Mathnium will read and execute the subsequent commands from the file named cfile.m. Moreover, if the name of a file to be executed has the suffix .m, you just have to type in its name without double quotes. Thus, if myfile is a currently undefined variable, the command

myfile 

will result in the execution the statements in the file myfile.m .

COMMENTS

All the input in a line following one of the three strings #, % or // is treated as comment. Multi-line comments commencing with the string /* are terminated by the string */.
>>a=2 % comment
>>b=10 # comment
>>c=a+b // comment
>>/*
>>Multi line comment
>>*/
>>c
   12

TERMINATING EXECUTION

A Mathnium session may be terminated by using the command exit or quit.

NUMBERS

In addition to integers and floating point numbers, Mathnium allows you to deal with complex numbers in a natural way. A complex number is specified by either appending the string _j to the imaginary part of the number, or by appending the imaginary part of the number to the string j_. Initially, the value of the variable i is set to be sqrt(-1), and so pre or post-multiplication by i may also be used to make any real valued expression imaginary. All the arithmetic operators valid for real numbers can be used with complex numbers as well.
>>2+j_3
               2 +               3i
>>2+j_3+6
               8 +               3i
>>(2+j_3)*(3+j_3)
              -3 +              15i
>>(5+j_2)/(3-j_2)
          0.8462 +          1.2308i

VARIABLE NAMES

There is no explicit limit on the number of chracters in identifiers (variable names) whose first character must be a letter which may be followed by numeric digits, letters, or underscors. The names may not begin with j_ or i_, as any name of the form j_xxxxx or i_xxxxx is interpreted as sqrt(-1)*xxxxx. Likewise, the variable names may not terminate with _j or _i, as any name of the form xxxxx_j or xxxx_i is interpreted as sqrt(-1)*xxxxx.

DEFINITION OF ARRAYS OF NUMBERS

Arrays are basic data types in Mathnium in that they can be used as operands for various ioperators. A one-dimensional array is just a row or a column of numbers. A two-dimensional array consists of rows of numbers, with each row consisting of the same number of elements. An array with p rows, each of which contains q numbers, is said to be a p by q array or an array containing p rows and q columns. Thus, a one-dimensional array is either a p by 1 or a 1 by q array. Two arrays are of the same size if they have the same number of columns and rows.

The simplest way to enter an array is to enclose it in square brackets, separating the rows either by a semi-colon or by a carriage return.

>>a=[2 3 4; 6 3 2]
>>a
          2          3          4
          6          3          2

>>b=[6 5 3 2]
>>b
          6          5          3          2

>>c=[9;7;6;3]
>>c
          9
          7
          6
          3

Multi-dimensional arrays may be defined by using the function cat which catenates a set of specified arrays along a particular dimension. If the first is input to the function is one, it stacks the successive elements in rows of a matrix, and if it is two, it stacks them in columns. Clearly, for higher values of the first input, the function cat stacks the elements along the specified dimension of a multi-dimensional array.
>>//catenate along the third dimension
>>cat(3,[2 3 ; 4 5],[7 8; 9 10])
[:,1:2, 1]
            2            3
            4            5
[:,1:2, 2]
            7            8
            9           10

>>//catenate along the fourth dimension
>>cat(4,[2 3 ; 4 5],[7 8; 9 10])
[:,1:2, 1, 1]
            2            3
            4            5
[:,1:2, 1, 2]
            7            8
            9           10

ARITHMETIC OPERATIONS WITH ARRAYS

The arithmetic operators

+  -  *  /  ^ 

can be used with the scalars (i.e., arrays with one row and one column) in the usual way for addition, subtraction, multiplication, division and exponentiation. However, you have to be a little bit more careful when using these operators with arrays.

A scalar can be added to an array of numbers by using the addition operator +. The result is an array obtained by adding the scalar to each element of the array used as the operand. The operators -, *, / and ^ can be used in a similar fashion with a scalar and an array as the operands.

The operators +, - and ^ can be used with two arrays of the same size as the operands. The result is obtained by the element-wise application of the operator to the corresponding elements of the two arrays, and is an array of the same size as the operands.

For the similar element-by-element multiplication or division with two arrays of numbers, you can use the operators .* and ./. The operator .^ corresponds to element-by-element exponentiation.
>>a=[2 3 4]
>>a+2
          4          5          6

>>a-2
          0          1          2

>>b=[6 7 8]
>>a+b
          8         10         12

>>a-b
         -4         -4         -4

>>a*2
          4          6          8

>>a.*b
         12         21         32

>>a./b // Right element by element division
       0.3333       0.4286          0.5

>>a.\b // Left element by element division
            3       2.3333            2

MATRIX MULTIPLICATION AND DIVISION

If the operands are arrays, the operators *, / and \ correspond to matrix multiplication, right division, and left division respectively. Thus, these operators can be used with arrays only if the number of rows in the second operand is the same as the number of columns in the first operand. The number of rows in the resulting array is the same as the number of rows in the first operand, and the number of columns in the resulting array is the same as the number of columns in the second operand.

Thus, a system of linear equations with a square coefficient matrix can be solved by using the matrix division operator. (If the divisor is not a square matrix, the least square solution of a system of equations is returned.)
>>a=[2 3 ; 8 9]
>>b=[7 9 10; 88 32 67]
>>a*b
        278        114        221
        848        360        683

>>a\b
         33.5          2.5         18.5
          -20       1.3333           -9

>>c=[6 2; 9 10;12 18]
>>c/a
      -6.3333       2.3333
      -0.1667       1.1667
            6            0

MATRIX TRANSPOSE

The operators .' and ' can be used for obtaining, respectively, the transpose and the Hermitian conjugate of matrices.
>>a=rand(2,3) // rand(2,3) returns a 2X3 matrix of random numbers.
>>a
       0.9897       0.1033       0.8171
       0.2327       0.0397       0.7688

>>a'
       0.9897       0.2327
       0.1033       0.0397
       0.8171       0.7688

>>a=a+j_rand(2,3) // the prefix j_  multiplies the expression by sqrt(-1).
>>a
       0.9897 +     0.4806i          0.1033 +     0.7386i          0.8171 +
0.3632i
       0.2327 +     0.6523i          0.0397 +     0.9303i          0.7688 +
0.1731i

>>a'
       0.9897 -     0.4806i          0.2327 -     0.6523i
       0.1033 -     0.7386i          0.0397 -     0.9303i
       0.8171 -     0.3632i          0.7688 -     0.1731i

>>a.'
       0.9897 +     0.4806i          0.2327 +     0.6523i
       0.1033 +     0.7386i          0.0397 +     0.9303i
       0.8171 +     0.3632i          0.7688 +     0.1731i

RELATIONAL AND LOGICAL OPERATORS

The following relational and logical operators are defined in Mathnium:
==  equal 
!= not equal 
> greater than 
>= greater than or equal 
< less than 
<== less than or equal 
~ (or !) unary negation 
|| (or |) logical or 
&& (or &) logical and 

The relational and logical operators accept following types of arguments:
(i) Scalars, 
(ii) An array and a scalar, 
(iii) Two arrays of the same size. 

If the operands for relational operators are scalars, they return 1 (for TRUE) if the relationship is satisfied by the operands; otherwise the result of these operations is zero (for FALSE). In other cases, the relational operators yield an array of zeros or ones obtained by pairwise comparison of each element of the operands. Note that if one of the arguments is a scalar and the other an array, the same scalar is compared with each of the elements of the array.

The logical operators behave in the same way as the relational operators except that they accept only integers and integer arrays as operands. The logical operators assume any nonzero integer to be TRUE and zero to be FALSE, and return zeros or ones.
>>[2 10 5]==2
          1          0          0

>>[3 8 ; 4 10] != [ 4 8 ; 3 10]
          1          0
          1          0

>>[2 7 0]&&0
          0          0          0

>>[2 7 0]&&1
          1          1          0

>>~[2 3 2]
          0          0          0

>>~[2 0 8]
          0          1          0

MATHNIUM FUNCTIONS

The Mathnium interpreter comes with close to four hundred functions for numerical computations, graphics, and input-output. The detailed descriptions of the functions are available in elsewhere. Here is a brief overview of the kinds of computations that that can be performed with Mathnium.

The functions for numerical computations cover the following areas:
What follow are some examples of the use of some of the functions for numerical computations.

LINEAR ALGEBRA

The first function that you should learn to use is rand. It allows you to conveniently generate arrays for input to various functions with which you want to play around in order to familiarize yourself with Mathnium.

The call rand(p,q) generates a p by q array of real numbers which are uniformly distributed between 0 and 1.
>>a=rand(4,4)
>>a
       0.5848       0.7915        0.296       0.4346
       0.1586       0.4306       0.5582       0.2974
       0.0235       0.3031       0.3388       0.5262
       0.7698       0.5108        0.039       0.6881

>>[c,d]=eig(a)
The last call returns the array of eigenvalues of the matrix a in the diagonals of the output d and the matrix of eigenvectors in the variable c. The call illustrates two basic features of the syntax for function calls in Mathnium: (i) the output arguments are separated from the inputs and (ii) functions can return multiple outputs. These features lead to considerable clarity in function definitions as you will begin to appreciate as you become more familiar with the system.

Let us try to see the accuracy of the eigenvalue computations by calculating the norm of a matrix whose elements should all be zero in an infinite precision calculation.
>>norm(a*c-c*d)
       1.3328e-015
The function svd computes the singular value decomposition of a matrix:
>>a=rand(100,100)		# Define a matrix 
>>[u,s,v]=svd(a)                # Compute the SD 
>>norm(u'*a*v-s) 		# Check the accuracy of the results  
       4.8356e-013

The function qr is one more in a long list of Mathnium functions for matrix decompositions. If called with three output arguments, it returns a lower triangular matrix r, a unitary matrix q and a permutation matrix e such that

q*r = a*e' 

Here is an example for the use of qr:
>>a=rand(4,4)
>>(q,r,e)=qr(a)
>>q'*q            // q is unitary.
            1 -2.7756e-017  1.1102e-016 -5.5511e-017
 -2.7756e-017            1 -4.1633e-017 -1.5613e-017
  1.1102e-016 -4.1633e-017            1 -8.3267e-017
 -5.5511e-017 -1.5613e-017 -8.3267e-017            1

>>r               // r is upper triangular.
      -1.2933      -0.5724      -0.8276      -0.1615
            0      -0.8898       0.0062      -0.1368
            0            0      -0.4123      -0.1434
            0            0            0       0.1931

>>q*r-a*e'        // This should vanish.
Scaled by 10^-16
       2.2204       2.2204       2.2204       0.5551
       0.2776      -1.1102       0.2776      -0.5551
       1.1102       0.5551       2.2204       0.5551
       2.2204            0       1.1102       0.1388
We give just one more example from linear algebra. The function simplex is an implementation of the simplex algorithm for solution of linear programming problems. The simplest of such problems are of the form
Given a row vector c, a column vector b 
and a matrix a, find the vector x that  
minimizes c*x subject to the constraint  a*x = b. 

Of course, simplex solves the more general problems also, but the following example is for the restricted problem:
>>a=[1 6 -1 0; 0 -3 4 1]
>>b=[2 ; 8]
>>c=[0 -2 4 0]
>>(x,cx)=simplex(a,b,c)
>>x                      // The solution
            0
       0.3333
            0
            9

>>cx                    // The optimal cost function
         -0.6667
>>a*x                   // This should be the same as b.
            2
            8

MATHEMATICAL FUNCTIONS

Most of the common mathematical functions are available in Mathnium. Many Mathnium functions have been defined such that they accept an array of numbers as input, and return a table of the values of the mathematical function for each element of the input array.
>>sin([0 pi/10 pi/4])
            0        0.309       0.7071

>>bessj(0,[0:.1:1])
Columns 1 through 8
            1       0.9975         0.99       0.9776       0.9604       0.9385
      0.912       0.8812
Columns 9 through 11
       0.8463       0.8075       0.7652

Note here that in the call to bessj, for calculating the Bessel function of the first kind of the order zero, the input [0:.1:1] corresponds to a table of numbers between zero and one, with increments of 0.1. This construct with colons can be used to conveniently generate an array with regularly increasing or decreasing values.

The input to many Mathnium functions may be in the range of inputs for which they return complex values. In fact, all the Mathnium trigonometric and hyperbolic functions, including the inverse functions, can be called with any real or complex numbers as the input so long as the corresponding value of the function is bounded.
>>x=[0:.2:2]
>>y=asin(x) % Some of the elements of the input are > 1.
>>y         %
Columns 1 through 3
            0                        0.2014                        0.4115

Columns 4 through 6
       0.6435                        0.9273                        1.5708

Columns 7 through 9
       1.5708 +     0.6224i          1.5708 +      0.867i          1.5708 +
 1.047i
Columns 10 through 11
       1.5708 +     1.1929i          1.5708 +      1.317i

>>sin(y)          % should be same as x
Columns 1 through 3
            0                           0.2                           0.4

Columns 4 through 6
          0.6                           0.8                             1

Columns 7 through 9
          1.2 +4.0617e-017i             1.4 +5.9995e-017i             1.6 +7.647
9e-017i
Columns 10 through 11
          1.8 +9.1644e-017i               2 +1.0606e-016i

>>norm(sin(y)-x)
       1.9145e-015

CALCULUS

The function series returns the Taylor series expansion of an expression in a single variable. The first input to this function is required to be the function (of a single variable) whose series is needed. Such arguments may be specified in in a number of ways including the following:
>>series(@(x)sin(x)/x,0)
            1            0      -0.1667            0

>>series(@(x)sin(x)/x,0,10) // Use 10 terms in the series to compute the expression.
Columns 1 through 8
            1            0      -0.1667            0       0.0083            0
    -0.0002            0
Column 9
  2.7557e-006
The first input to series is function object, defined here as an inline function. The second input is the value of the variable about which you want to expand the function. In this case the function series yields the limit of sin(x)/x as x goes to zero.

The inputs to the function deriv are also similar. The function calculates the first derivative of a scalar function in a single variable.


>>deriv(@(x)(tan(x^2+2*x)+cos(x)),pi/4)
          9.9639

Let us now check the result of deriv .
>>x=pi/4
>>b=cos(x^2+2*x)
>>(2*x+2)/(b*b)-sin(x)
          9.9639

NONLINEAR EQUATIONS AND FUNCTION MINIMIZATION

The function fzero yields the root of an expression in single variable:
>>y=fzero(@(x)tan(x)*x-1,rand(1,1))
>>y
          0.8603
>>y*tan(y)-1
      -3.2479e-008
The last input is our guess of the root. Here is a check of the accuracy of the result:
>>y*tan(y)-1
      -3.2479e-008
As with many other functions, you can use optional input arguments in fzero. The last of the three optional inputs allowed in fzero controls the criterion for convergence of the root. Let us use our own value for this input.
>>y=fzero(@(x)tan(x)*x-1,rand(1,1),null,null,1.e-12)
>>y*tan(y)-1
      -4.2612e-010
You can see that the accuracy of the root has improved. You may also have noted that since we did not want to use the first two of the three allowed optional inputs to froot, we just used nulls in their places.

Here is one more example of the use of froot.
>>y=fzero(@(x)bessj(0,x)-bessj(1,x),rand(1,1))
>>bessj(0,y)-bessj(1,y)
         0
The function fminbnd minimizes a function of a single variable. The first input to fminbnd is the function to be minimized, and the two subsequent inputs specify the limits of The third input is the name of the function which we want to minimize:
>>(f,v)=fminbnd(@(x)bessj(0,x),0,10)
>>f   // The estimated minimum value
         -0.4028
>>v   // The point at which the minimum is estimated to lie.
          3.8317
We can now check the value of the function to be minimized at points near the purported minimum:
>>bessj(0,[.99 1 1.01]*v)
      -0.4025      -0.4028      -0.4025
You can solve a set of equations by the Newton's method by using the function newton, and use fminbfgs to minimize a function of several variables. Both newton and fminbfgs require you to define a function that also returns the Jacobian/gradient of the appropriate quantities. In the following example, the simple equations
x[1]*x[2] = 10 
x[1]+x[2] = 7 

are solved by using newton.
>>function [y,dy]=equations(x)
>      y=[x[1]*x[2]-10;x[1]+x[2]-7]
>      dy=[x[2] x[1];1 1]
>end
>>x=newton(@equations,rand(2,1))
>>x
          5
          2
If the Jacobian of the set of nonlinear equations to be solved is not available, the function broyden may be used.
>>//Solve the two nonlinear equations by the broyden's method.
>>broyden(@(x)[x[1]*x[2]-10;x[1]+x[2]-7],rand(2,1))
>>x
          5
          2

For more details, see the descriptions of these functions.

FAST FOURIER TRANSFORMS

The radix-2 Fast-Fourier transform algorithm is implemented in the function fft. You can calculate the inverse transform by using ifft.
>>x=[0:.1:1024]';
>>//Construct a time series with three dominant frequencies 
>>//added to a sequence of random numbers.
>>y=sin(6pi*x)+3*cos(7.6pi*x)+2sin(9.8pi*x)+randn(size(x))*.1
>>//Plot the time series.
>>subplot(2,1,1)
>>plot(x,y)
>>title("Time Series")
>>//Compute the transfrom.
>>z=fft(y)
>>n=size(z,1);
>>//Plot the amplitudes.
>>subplot(2,1,2)
>>plot(abs(z[1:n/2+1]))
>>title("FFT")


The functions fft2d and ifft2d calculate the Fast-Fourier transform pairs for two-dimensional arrays.

ORDINARY DIFFERENTIAL EQUATIONS

The function rungekutta is an implementation of the adaptive Fourth-order method for the solution of ordinary differential equations.

In the example shown here, solution is obtained for the coupled equations governing the Bessel functions of the first kind of order zero and one at twenty points starting at 1, with the distances between the points being 0.1.
>>r0=1                              # The initial value of the independent coordinate
>>y0=[bessj(0,r0);bessj(1,r0)]      # The initial condition
>>h=.1;                             # The distance between the sampling points
>>(r,y)=rungekutta(@(r,y)[-y[2];y[1]-y[2]/r],r0,y0,h,20) # Solve
>>yb=[bessj(0,r) bessj(1,r)]        # Equivalent values from library functions
>>norm(y-yb)                        # difference
       3.6460e-007

GRAPHICS

As exemplified above for the function fft, Mathnium provides a comprehensive library for graphics. The simplest function to use is plot, which draws an x-y plot of the inputs. As shown below, the plots can be easily annotated and saved as images.
>>x=[0:.1:10]
>>y=bessj(0,x);
>>plot(x,y);
>>xlabel(" X ")
>>ylabel(" Y ")
>>title({"Bessel function","Order 0"});
>>// Save as a 700 X 550 image in a gif file.
>>saveFigure("bessj0.gif",700,550);


Functions are available for drawing the following types of graphs:

PROGRAMMING IN MATHNIUM

Obviously, the usefulness of Mathnium as a tool for numerical computing can be readily extended by using it as a programming language. The usual control structures such as while and for loops, as well as the if-block and switch block are all supported. Further, functions as well as classes can be defined to implement algorithms not already available in the library bundled with Mathnium.

Here is the function for computing the Fibonocci sequence. It does not use recursive calls, and stores the results in a static array for quick retrieval of the result in case it has already been computed.
function fib(n)
   if(n<=2);return n;end 
   // Static store for the elements of the sequence that have been already computed.
   static scratch
   // initialization
   if(scratch==null)
     scratch[1]=1
     scratch[2]=2
   end
   // number of elements already available
   m=numel(scratch)
   // extend the static store if needed
   for i=m+1:n
     scratch[i]=scratch[i-1]+scratch[i-2]
   end
   // return the element corresponding to the input argument.
   return scratch[n]
end
The following fractal tree was drawn using the call fractaltree(6). The complete function is shown below the image.

>// Display the file containing the function fractaltree
>>showfile("fractaltree.m")
/*
        Draws a fractal tree of depth n.
        Usage:
                fractaltree(n)

        n should be less than 6 as otherwise it may take too long.
*/
function fractaltree(n)
   branch([0,0],[0 .5],1,pi/2,0,n)
   hold off
end
function branch(x,y,len,theta,n,nmax)
/*
        Draws a line between the points (x[1],y[1]) and (x[2],y[2]).
        and five branches of length len/2. One branch is an extension
        of the line itself and the other four branches are at -60,-30,30,
        and 60 degrees with respect to the main branch
        starting from the point (x[2],y[2]).

        len: twice the  length of the branches
        theta: the angle between the line and the x-axis.
        n:  recursion depth
        nmax: maximum recursion depth
*/
        // Maximum depth of recursion
        if(n>=nmax);return;end
        // colors
        static red=java.awt.Color.red;
        static green=java.awt.Color.green;
        // clear the current graphics window
        if(n==0); cleargraph();end
        // choose red or green color
        color=red;if(n>nmax/2);color=green;end
        // The main branch
        line([x[1] y[1] x[2] y[2]],'color',color);
        // Initialization
        if(n==0)
                axis([-1,1,0,2]); % Scaled limits of the drawing area.
                axis('square');   // Aspect ratio of the drawing area.
                notics();        # Suppress ticmarks and tic labels.
                hold on;          // Do not erase
        end
        // Now draw the five branches
        len/=2;
        t1=theta-pi/2;
        for i=1:5
                t1=t1+pi/6
                dx=len*cos(t1);
                dy=len*sin(t1);
                branch([x[2],x[2]+dx],[y[2],y[2]+dy],len,t1,n+1,nmax);
        end
end
>>// Draw a fractaltree of depth 6 using the function shown above.
>>fractaltree(6)
>>// Save the figure as an image file.
>>saveFigure("fractaltree.gif",800,600)