Monday, December 6, 2010

Ito Calculus with Tensorial

Here is a short example on how to use the Tensorial package for Ito calculus (stochastic calculus)

ito.pdf

Numerical Linear Algebra and Quantum Control

Here is a poster about some of my recently published research

Poster pdf

Sunday, December 5, 2010

MathLink example 4: Trace of a complex matrix

In this example, I introduce the ability to treat complex numbers.

-------------------------------------------------------
myTrace.c
-------------------------------------------------------

#include<stdio.h>
#include<string.h>

#include "mathlink.h"
extern void myTrace( void );

void myTrace( void )
{

 int i,j,n;
  
 double *matrix;

 int *dims;
 char **heads;
 int rank; 

 if( !MLGetReal64Array(stdlink,&matrix,&dims,&heads,&rank) ){ 
 // unable to read data

  printf(" MLGetReal64Array error reading data \n");
 };

 n = dims[0];

 double resum=0.;
 double imsum=0.;

 double sum[2];

 int outdims[1];

 char **outheads;
 int outrank = 1; 

 
 for(i=0;i<n;i++){

 resum = resum + matrix[2*i*n + 2*i];

 imsum = imsum + matrix[2*i*n + 2*i + 1];
 }

 
 MLReleaseReal64Array(stdlink,matrix,dims,heads,rank);
 
 MLPutFunction(stdlink,"Complex",2);

 MLPutFloat(stdlink,resum);
 MLPutFloat(stdlink,imsum);


}


#if WINDOWS_MATHLINK

#if __BORLANDC__
#pragma argsused
#endif

int PASCAL WinMain( HINSTANCE hinstCurrent, HINSTANCE hinstPrevious, LPSTR lpszCmdLine, int nCmdShow)
{

 char  buff[512];
 char FAR * buff_start = buff;

 char FAR * argv[32];
 char FAR * FAR * argv_end = argv + 32;

 hinstPrevious = hinstPrevious; /* suppress warning */

 if( !MLInitializeIcon( hinstCurrent, nCmdShow)) return 1;

 MLScanString( argv, &argv_end, &lpszCmdLine, &buff_start);
 return MLMain( (int)(argv_end - argv), argv);
}

#else

int main(int argc, char* argv[])
{

 return MLMain(argc, argv);
}

#endif

-----------------------------------------------------
myTrace.tm
-----------------------------------------------------
:Begin:
:Function:       myTrace

:Pattern:        myTrace[ L_List]
:Arguments:      { L  }
:ArgumentTypes:  { Manual  }
:ReturnType:     Manual

:End:

:Evaluate: myTrace::usage = "myTrace[M] gives the trace of complex matrix M"

------------------------------------------------------
Makefile
------------------------------------------------------
MPREP = /usr/bin/mprep
CXX = /usr/bin/c++

myTrace : myTracetm.c myTrace.c
 ${CXX} myTracetm.c myTrace.c -o myTrace -lML64i3 -lm -lpthread -lrt -lstdc++

myTracetm.c : myTrace.tm
 ${MPREP} myTrace.tm -o $@


---------------------------------------------------
myTrace.nb
---------------------------------------------------

link = Install["./myTrace"]

(m = {{1., 2., 1.}, {2., 2., 2.}, {3., 3., 3.}} + I) // MatrixForm

myTrace[m]

Saturday, December 4, 2010

MathLink example 3: Trace of a real matrix

In this example we calculate the trace of a matrix of floating point numbers (double type in C)

---------------------------------------------------
myTrace.c
---------------------------------------------------

#include "mathlink.h"
extern double myTrace( void );

double myTrace( void )
{
 int i,j,n;
  
 double *matrix;
 int *dims;
 char **heads;

 int rank; 

 if( !MLGetReal64Array(stdlink,&matrix,&dims,&heads,&rank) ){ 
 // unable to read data

  return 0.;
 };

 n = dims[0];

 double sum=0;
 for(i=0;i<n;i++) sum = sum + matrix[i*n+i];
 
 MLReleaseReal64Array(stdlink,matrix,dims,heads,rank);

 return sum;
}

#if WINDOWS_MATHLINK

#if __BORLANDC__
#pragma argsused
#endif

int PASCAL WinMain( HINSTANCE hinstCurrent, HINSTANCE hinstPrevious, LPSTR lpszCmdLine, int nCmdShow)
{

 char  buff[512];
 char FAR * buff_start = buff;

 char FAR * argv[32];
 char FAR * FAR * argv_end = argv + 32;

 hinstPrevious = hinstPrevious; /* suppress warning */

 if( !MLInitializeIcon( hinstCurrent, nCmdShow)) return 1;

 MLScanString( argv, &argv_end, &lpszCmdLine, &buff_start);
 return MLMain( (int)(argv_end - argv), argv);
}

#else

int main(int argc, char* argv[])
{

 return MLMain(argc, argv);
}

#endif

----------------------------------------------
myTrace.tm
----------------------------------------------

:Begin:
:Function:       myTrace

:Pattern:        myTrace[ L_List]
:Arguments:      { L  }
:ArgumentTypes:  { Manual  }
:ReturnType:     Real

:End:

:Evaluate: myTrace::usage = "myTrace[M] gives the trace of matrix M"

---------------------------------------------
Makefile
---------------------------------------------

MPREP = /usr/bin/mprep
CXX = /usr/bin/c++

myTrace : myTracetm.c myTrace.c
 ${CXX} myTracetm.c myTrace.c -o myTrace -lML64i3 -lm -lpthread -lrt -lstdc++

myTracetm.c : myTrace.tm
 ${MPREP} myTrace.tm -o $@

Wednesday, November 24, 2010

Catenary

The catenary is the curve that is naturally formed when a chain of constant linear density hangs from two fixed points. This is a typical example on variational calculus showing the use of Lagrange multipliers. Next, I develop the problem entirely in Mathematica

Catenary
Catenary(pdf)

You may need MathematicaPlayer if you do not have Mathematica.

This problem can be generalized to find the surface with minimum energy, which is shown next


MathLink Example 2: Sum a list of floating point numbers

This example considers a sum of of floating point numbers.
----------------------------------------------------------------
addRealList.c
----------------------------------------------------------------

#include "mathlink.h"
extern double addRealList( void );



double addRealList( void )
{
 int i,n;

 double *list; 

 MLGetReal64List(stdlink,&list,&n);

 double sum=0;
 for(i=0;i<n;i++) sum = sum + list[i];

 
  MLReleaseReal64List(stdlink,list,n);

 return sum;
}



#if WINDOWS_MATHLINK

#if __BORLANDC__
#pragma argsused
#endif

int PASCAL WinMain( HINSTANCE hinstCurrent, HINSTANCE hinstPrevious, LPSTR lpszCmdLine, int nCmdShow)
{

 char  buff[512];
 char FAR * buff_start = buff;

 char FAR * argv[32];
 char FAR * FAR * argv_end = argv + 32;

 hinstPrevious = hinstPrevious; /* suppress warning */

 if( !MLInitializeIcon( hinstCurrent, nCmdShow)) return 1;

 MLScanString( argv, &argv_end, &lpszCmdLine, &buff_start);
 return MLMain( (int)(argv_end - argv), argv);
}

#else

int main(int argc, char* argv[])
{

 return MLMain(argc, argv);
}

#endif


------------------------------------------------------
addRealList.tm
------------------------------------------------------

:Begin:
:Function:       addRealList

:Pattern:        addRealList[ L:{___Real}]
:Arguments:      { L  }
:ArgumentTypes:  { Manual  }
:ReturnType:     Real

:End:

:Evaluate: addRealList::usage = "addRealList[x] gives the sum of the elements on x, given that they are real (double floating numbers)"

------------------------------------------
Makefile
------------------------------------------
MPREP = /usr/bin/mprep
CXX = /usr/bin/c++

addRealList : addRealListtm.c
 ${CXX} addRealListtm.c addRealList.c -o addRealList -lML64i3 -lm -lpthread -lrt -lstdc++

addRealListtm.c : addRealList.tm
 ${MPREP} addRealList.tm -o $@

clean:

 rm addRealListtm.c
 rm addRealList
 

--------------------------------------------------------
addRealList.nb
--------------------------------------------------------

link = Install["./addRealList"]

?addRealList

addRealList[{1., 5., 6.}]

Tuesday, November 23, 2010

MathLink Example 1: Sum a List of Integers

This example passes a list of integers, performs the sum and returns an integer.
We must note that the length of the list is not type int but type long.

-------------------------------------------------------------------------
addList.c
-------------------------------------------------------------------------

#include "mathlink.h"
extern double addList( double *list, long n);



double addList( double *list, long n)
{

 int i;
 double sum=0;
 for(i=0;i<n;i++) sum = sum + list[i];

 
 return sum;
}



#if WINDOWS_MATHLINK

#if __BORLANDC__
#pragma argsused
#endif

int PASCAL WinMain( HINSTANCE hinstCurrent, HINSTANCE hinstPrevious, LPSTR lpszCmdLine, int nCmdShow)
{

 char  buff[512];
 char FAR * buff_start = buff;

 char FAR * argv[32];
 char FAR * FAR * argv_end = argv + 32;

 hinstPrevious = hinstPrevious; /* suppress warning */

 if( !MLInitializeIcon( hinstCurrent, nCmdShow)) return 1;

 MLScanString( argv, &argv_end, &lpszCmdLine, &buff_start);
 return MLMain( (int)(argv_end - argv), argv);
}

#else

int main(int argc, char* argv[])
{

 return MLMain(argc, argv);
}

#endif


---------------------------------------------------------------------
 addList.tm
---------------------------------------------------------------------
:Begin:
:Function:       addList

:Pattern:        addList[ L_List]
:Arguments:      { L  }
:ArgumentTypes:  { RealList  }
:ReturnType:     Real

:End:

:Evaluate: addList::usage = "addList[x] gives the sum of the elements of the list x"

---------------------------------------------------------------------
Makefile
---------------------------------------------------------------------
MPREP = /usr/bin/mprep
CXX = /usr/bin/c++

BINARIES = addList

all : $(BINARIES)

addList : addListtm.c addList.c
 ${CXX}  addList.c addListtm.c  -lML64i3 -lm -lpthread -lrt -lstdc++ -o $@

addListtm.c : addList.tm
 ${MPREP} addList.tm -o addListtm.c

-------------------------------------------
addList.nb
-------------------------------------------

link = Install["./addList"]

addList[{1, 5, 6}]