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:
- 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
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:
- 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-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:
- 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.
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)