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.

>>2+3, 3^2 # A simple statement 5 9 >>s=0 >>for i=1:10 # A control structure > s=s+i > end >>s 55It 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 myfilewill result in the execution the statements in the file

```
myfile.m .
```

`#`

,
`%`

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

`exit`

or
`quit`

.
```
_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

`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`

.
```
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 3Multi-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

+ - * / ^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

```
*, /
```

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

```
.'
```

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

== equal != not equal > greater than >= greater than or equal < less than <== less than or equal ~ (or !) unary negation || (or |) logical or && (or &) logical andThe 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

The functions for numerical computations cover the following areas:

- Linear Algebra
- Integration, Differentiation and Taylor Series Expansion
- Evaluation of Mathematical Functions
- Solution of Nonlinear Equations and Function Minimization
- Ordinary Differential Equations
- One and Two-Dimensional Fast Fourier Transforms
- Data Analysis: Sorting and Simple Statistics
- Polynomial Arithmetic

```
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-015The 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-013The 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.1388We 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

>>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.7652Note 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

```
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:
- Use the construct
`$func`

or`@func`

to sepecify an existing function`func`

as the argument. - Use an
*inline*function of the form`@(x)expr`

where`expr`

is an expression defined for the argument`x`

as, for example, in`@(x)sin(x)*x-1`

.

>>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-006The 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.9639Let 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

```
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-008The last input is our guess of the root. Here is a check of the accuracy of the result:

>>y*tan(y)-1 -3.2479e-008As 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-010You 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) 0The 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.8317We 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.4025You 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] = 7are solved by using

```
newton.
```

If the Jacobian of the set of nonlinear equations to be solved is not available, the function>>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)) >>x5 2

`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)) >>x5 2

For more details, see the descriptions of these functions.

```
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.
```
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

`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:

- x-y line and scatter plots.
- Two-D and three-d bar charts.
- Two-D and three-d pie charts.
- Contour plots.
- Mesh plots and surface plots of functuons dependent on two variables.

```
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] endThe following

`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)