InterCall

A Mathematica package that allows interactive calling of external code

Linking Mathematica with the IMSL and NAG libraries

InterCall is a Mathematica package that provides:

With InterCall you can:

For further information contact: intercall_forward@wri.com or

Analytica                                                              Dr. Terry Robb
    PO Box 343                                                            Mathematics Department
    Subiaco 6008                                                          Monash University
    Perth WA                                                               Clayton 3168
    Australia                                                                Melbourne VIC
    Fax/Phone +61 9 388 2590                                       Australia
    Email: analytic@earwax.pd.uwa.oz.au                Email: tdr@vaxc.cc.monash.edu.au

Table of Contents

Initialisation

Quadrature

Special functions

Interpolation

Root finding

Differential equations

Linear algebra

Eigensystems

Sparse matrices

Some NAG Examples

Using InterCall to call the NAG library

This notebook contains a sample Mathematica session in which InterCall is used to access various routines in the NAG subroutine library. The examples are intended to indicate how easily InterCall integrates the symbolic capabilities of Mathematica with the numeric features of an external library such as NAG.

The main InterCall command illustrated in these examples is GetDefault although there are other InterCall commands called SetDefault and AddDefault. The find command is a simple searching utility that can be used to help select an appropriate IMSL or NAG routine.

Initialisation

Load the InterCall package

<<InterCall.m

Connect to a remote computer that has the NAG subroutine library

InterCall["vaxc.cc.monash.edu.au -l tdr"];

Import various NAG routines from that remote computer

<<InterData.m;
Import[{D01AJF, D01AKF, D01AMF, S17AGF, E01BAF, E02BBF,
        F01BTF, F04AYF, F04EAF, F01AXF, C05ADF, C05NBF,
        X02AAF, F02AGF, F02ABF, D02BAF, D02EBF, D02BBF,
        D02GAF, F02BBF, F01MCF, F04MCF, F01BRF, F04AXF}]

Quadrature

Integrate the oscillatory function x cos(30x) sin(x) from 0 to 29

find["quadrature" && "finite" && "one" && "oscil" && "(NAG)"]

GetDefault[d01akf]

Compute the integral, and save the sampled points in a list

samples = {};
d01akf[ Function[x, AppendTo[samples,x];
                    x Cos[30x] Sin[x]], 0, 2Pi ]

Find the absolute error in the last integral

d01akf`abserr

Plot the points where the NAG routine d01akf has sampled the function

ListPlot[samples]; samples=.

Plot the oscillatory function

Plot[ x Cos[30x] Sin[x], {x,0,2Pi}, PlotPoints->60 ];

Integrate 1/((x+1)Cx) over the semi-infinite region 0 to 0

find[ "quadrature" && "infinite" && "one" &&
       Not["singular"] && "(NAG)" ]

GetDefault[d01amf]

Compute the integral (here +1 means positive infinity)

d01amf[ Function[x,1/((x+1)Sqrt[x])], 0, +1, $EPSREL->10^-7 ]

Check the answer using Mathematica's  NIntegrate function

% - NIntegrate[ 1/((x+1)Sqrt[x]), {x,0,+Infinity} ]

Integrate a badly-behaved function

find["quadrature" && "finite" && "one" && "bad" && "(NAG)"]

Plot the badly-behaved function

Plot[ x*Sin[30x]/Sqrt[1-(x/(2Pi))^2], {x,0,2Pi},
      PlotPoints->60 ];

Compute the integral

d01ajf[ Function[x, x*Sin[30x]/Sqrt[1-(x/(2Pi))^2]], 0, 2Pi,
        $EPSREL->10^-9 ]

Find the absolute error in the last integral

d01ajf`abserr

Special functions

Call the Airy function Ai(x), from the NAG library

find["Airy" && "(NAG)"]

GetDefault[s17agf]

Plot a graph of the function

Plot[ s17agf[x], {x,-3,3} ];

Interpolation

Fit a cubic spline through a set of points

find["cubic" && "spline" && "interpo" && "(NAG)"]

GetDefault[e01baf]

find["cubic" && "spline" && "fitted" && "E02BAF"]

GetDefault[e02bbf]

Compute the spline coefficients of the Airy function evaluated at some nodes

nodes = {-3., -2., -1.5, -1., -0.5, 0., 1., 2., 3.};
fvals = AiryAi[nodes];
e01baf[nodes,fvals];

Plot a graph of the fitted cubic spline

Plot[ e02bbf[x], {x,Min[nodes],Max[nodes]},
      Epilog->{PointSize[0.02], Point /@
               Transpose[{nodes,fvals}]} ]; nodes=.; fvals=.

Root finding

Solve a non-linear equation

find["solution" && "non-linear" && "system" && "(NAG)"]

GetDefault[c05nbf]

Consider the Boundary Value Problem: 3 x'' + 2 x' + 4 x2 = 2, with x=0 at the ends. By discretising the equation we get the non-linear system:

PxiP1 + (3P2xi)xi P 2xi+1 = P1,

with x0=0 and xn=0.

Here is a subroutine, required by c05nbf, that represents the system

myfcn = Function[{n,x,fvec,iflag}, Block[{xL,xR},
           xL = Prepend[ Drop[x,-1], 0 ]; (* x[i-1] *)
           xR =  Append[ Drop[x, 1], 0 ]; (* x[i+1] *)
           fvec = -xL + (3-2*x)*x - 2*xR + 1 ] ];

Solve the system using c05nbf

An initial guess is that x is constant and equal to P1. This we will approximate with 50 points on the interval [1,50].

root = c05nbf[ myfcn, Table[-1,{50}] ];
ListPlot[ Join[{0},root,{0}], PlotJoined->True,
          PlotRange->All ]; myfcn=.; root=.

Find a zero of a continuous function

find["zero" && "one" && "varia" && "(NAG)"]

GetDefault[c05adf]

Here is a function that has a zero on the interval [0,1]

func = Exp[#] - 3*# & ;
Plot[ func[x], {x,0,1} ];

Search the interval [0,1] to find a root, using the Bus & Dekker algorithm

c05adf[0,1,func]

(* func[%] //Chop *) func=.

Differential equations

Solve an initial value problem

find[ "system" && "ordi" && "diff" && "equ" && "initial" &&
       Not["stiff"] && "Runge" && "simple" && "(NAG)" ]

GetDefault[d02baf]

Consider the projectile equations:

dh/dx = tan(ϕ),
             dϕ/dx = Pg/v2,
             dv/dx = Pg tan(ϕ)/v P d sec(ϕ) v

where h is the height of the projectile, ϕ its angle above the horizon, v its tangential velocity, and x its distance down range.

Here is a subroutine, required by d02baf, that represents the system

mysub = Function[{x,y,f}, Block[{f1,f2,f3,h,phi,v,
            grav=0.032,drag=0.02},{h,phi,v} = y;
          f1 = Tan[phi];
          f2 = -grav/v^2;
          f3 = (f1*f2 - drag/Cos[phi])*v;
          f = {f1,f2,f3}; (*Print["x=",x];*) ] ];

Compute (h,ϕ,v) at range x=8, starting with h=0, ϕ=9/5, and v=0.5 at x=0

x = 0; y = {0,Pi/5,0.5}; xend = 8;
d02baf[x,xend,y,mysub,$TOL->0.0001]

Set a new default to make d02baf easier to call

d02baf[0,0,{0,Pi/5,0.5},mysub];
SetDefault[ d02baf, $X->$X, $Y->$Y:>Out, $FCN->$FCN ];

Compute a table of the solution from x=0 to x=8, and plot a graph of h(x)

The table can be computed much more efficiently using d02bbf. Unlike d02bbf, the routine d02baf has to re-determine an initial step size at each new call and this is expensive. The routine d02bbf can save up a list of the solution. See the example for d02ebf in the next section.

points = Table[ d02baf[x], {x,0,8,0.5} ];
ListPlot[ First /@ points, PlotJoined->True ]; points=.

SetDefault[d02baf]; mysub=.; x=.; y=.; xend=.

Solve an initial value problem that is stiff

GetDefault[d02ebf]

Consider the stiff system:

y1' = P 4%10P2 y1 + 104 y2 y3,
           y2' =  4%10P2 y1 P 104 y2 y3 P 3%107 y22,
           y3' = 3%107 y22

Here are two subroutines, required by d02ebf, that represent the system

The subroutines mysub and myder compute the RHS and its Jacobian for a given x and y.

mysub = Function[{x,y,f}, Block[{f1,f2,f3,y1,y2,y3},
         {y1,y2,y3} = y;
         f1 = -0.04 y1 + 10^4 y2 y3;
         f2 = -f1 - 3*10^7 y2^2;
         f3 = 3*10^7 y2^2;
         f = {f1,f2,f3}; (*Print["x=",x];*) ] ];
myder = Function[{x,y,pw}, Block[{y1,y2,y3}, {y1,y2,y3}=y;
         pw = {{-0.04, 10^4 y3,              10^4 y2},
               { 0.04,-10^4 y3 - 6*10^7 y2, -10^4 y2},
               { 0   ,           6*10^7 y2,  0      }}; ]];

Here is a subroutine, required by d02ebf, that outputs the solution

The subroutine myout decides at what x nodes the solution is to be computed, and we can also save the value of that solution in a list called soln.

myout = Function[{x,y}, AppendTo[soln, {x,y}];
         x = First[xnodes]; xnodes = Drop[xnodes,1]; ];

Compute the solution at some nodes, with initial condition y=(1,0,0) at x=0

In the following MPED = 0 means use the numerical Jacobian, and not myder. If MPED = 1, then d02ebf will use the analytic Jacobian defined in myder.

x0 = 0; y0 = {1,0,0}; xend = 10;
xnodes = {.0001,.001,0.01,.1,.3,1.,3.,5.,7.,10.,11.}; soln={};
d02ebf[x0,xend,y0,mysub,myder,myout,
       $MPED->0, $IRELAB->2, $TOL->0.0001];

Plot the solution

ListPlot[ {First[#],#[[2,1]]}& /@ soln, PlotJoined->True ];
ListPlot[ {First[#],#[[2,2]]}& /@ soln, PlotJoined->True ];
ListPlot[ {First[#],#[[2,3]]}& /@ soln, PlotJoined->True ];

mysub=.; myder=.; myout=.; soln=.; xnodes=.; x0=.; y0=.; xend=.

Solve a boundary value problem

find[ "system" && "ordi" && "diff" && "equ" &&
      "boundary" && "simple" && "(NAG)" ]

GetDefault[d02gaf]

Consider the boundary value problem:

y''' = P y y'' P ' (1 P y' 2),

y'(0) = y(0) = 0,
            y'(10) = 1

We wish to solve the above for '=0.2 on the interval [0,10] using at most 40 grid points. First we solve the system for '=0 with just 26 equally spaced grid points to find a suitable mesh, and then we solve again using that mesh and '=0.2. Read the entry for d02gaf in the NAG manual to find out the idea behind this example.

Here is a fortran subroutine, required by d02gaf, that represents the system

An auxillary subroutine called SETBETA is also written to allow Mathematica to modify the value of ' through a fortran common block.

Import the above two fortran subroutines into Mathematica

AddDefault[ MYSUB[$X,$Y,$F:>Out], S[R,R[3],R[3]], "mysub.o" ];
AddDefault[ SETBETA[$MYBETA], S[R], "mysub.o" ];
Import[{MYSUB,SETBETA}];

Solve the system for ' = 0.0, by guessing y(10)=10 and y''(0)=y''(10)=0.

The parameter u gives the boundary values for y, y' and y'' at the ends x=0 and x=10, and v specifies whether those values are exact or just guesses (0 for exact, 1 for a guess).

a = 0; b = 10;
u = {{0,b},{0,1},{0,0}}; v = {{0,1},{0,0},{1,1}};
n = 26; x = a + (b-a)/(n-1)*(Range[n]-1)//N;
setbeta[0.0];     d02gaf[u,v,a,b,MYSUB,x, $MNP->40];

Plot a graph of the solution y(x) with ' = 0.2 using the mesh X used by ' = 0.0

We use the spline routines e01baf and e02bbf to plot a smooth curve through the solution.

setbeta[0.2]; ans=d02gaf[u,v,a,b,MYSUB,0,
                        $X->$X, $NP->$NP, $MNP->40];
e01baf[ ans[[1]], First /@ ans[[2]] ];
Plot[ e02bbf[xx], {xx,a,b} ]; u=.;v=.;a=.;b=.;x=.;n=.;ans=.

Linear algebra

Solve the linear system Ax=b by LU-factorisation of the matrix A

find[" LU" && "factori" && "real" && Not["sparse"] && "(NAG)"]

GetDefault[f01btf]

find[ "solution" && "linear equa" && "factori" && "approx" &&
      "real" && Not["symmetric"] && Not["band"] &&
       Not["sparse"] && "(NAG)" ]

GetDefault[f04ayf]

LU-factorise a matrix A

a = Table[ Random[], {3}, {3} ];
f01btf[a];

Having factorised A, solve for some right hand side b,

b = Table[ Random[], {3}, {1} ];
f04ayf[b]

(* a.% - b //Chop *) a=.; b=.

Solve a tri-diagonal system

find["linear equa" && "tri" && "diagonal" && "(NAG)"]

GetDefault[f04eaf]

Here is a tri-diagonal matrix in Mathematica

a={{1,6,0,0,0},
   {6,2,5,0,0},
   {0,7,3,4,0},
   {0,0,8,4,3},
   {0,0,0,9,5}};

NAG represents the above tri-diagonal matrix by its leading diagonal {1,2,3,4,5}, its upper diagonal {_,6,5,4,3} and its lower diagonal {_,6,7,8,9}. The first element in the upper and lower diagonal is arbitrary as it is not used by f04eaf.

Solve the tri-diagonal system Ax=b

b={3,4,5,6,7};
f04eaf[ {1,2,3,4,5},{0,6,5,4,3},{0,6,7,8,9}, b ]

(* a.% - b //Chop *) a=.; b=.

QR-factorise a non-square matrix

find[" QR" && "factori" && "(NAG)"]

Eigensystems

Eigenvalues and eigenvectors

find[ "eigen" && "values" && "vectors" &&
      "real" && "symmetric" && "(NAG)" ]

GetDefault[f02abf]

Compute the eigenvalues and eigenvectors of a symmetric matrix

a=Table[ Switch[i-j,1,-1,0,2,-1,-1,_,0], {i,1,20},{j,1,20} ];
{eigvals,eigvecs} = f02abf[a];

Plot a graph of the first 3 eigenfunctions

plots = ListPlot[Join[{0},#,{0}], PlotJoined->True,
                   DisplayFunction:>Identity]& /@
        Take[eigvecs,3];
Show[ plots, DisplayFunction :> $DisplayFunction ];
eigvals=.; eigvecs=.; plots=.; a=.

Other examples

find[ "eigen" && "values" && "vectors" &&
       Not["symmetric"] && "real" && "(NAG)"]

Sparse matrices

Solve Ax=b for A a variable-bandwidth symmetric matrix

find["variable" && "band" && "(NAG)"]

GetDefault[ f01mcf ]

GetDefault[ f04mcf ]

Here is a band symmetric matrix in Mathematica, with a variable bandwidth

a={{ 1, 2, 0, 0, 5, 0},
   { 2, 5, 3, 0,14, 0},
   { 0, 3,13, 0,18, 0},
   { 0, 0, 0,16, 8,24},
   { 5,14,18, 8,55,17},
   { 0, 0, 0,24,17,77}};

Here is NAG's way of representing the above variable bandwidth matrix

aband={{1},{2,5},{3,13},{16},{5,14,18,8,55},{24,17,77}};

NAG represents a variable bandwidth symmetric matrix by its lower triangular bands.

Factorise the variable bandwidth symmetric matrix A

f01mcf[ Flatten[aband], Map[Length,aband] ];

Having factorised A, solve for some right hand side b,

b={{1},{2},{3},{4},{5},{6}};
f04mcf[b]

(* % - Inverse[N[a]].N[b] //Chop *) a=.; b=.; aband=.

Solve Ax=b for A a sparse matrix

find["sparse" && Not["symmetric"] && "real" && "(NAG)"]

GetDefault[ f01brf ]

GetDefault[ f04axf ]

Here is a sparse matrix (with lots of zero's) in Mathematica

a={{ 5, 0, 0, 0, 0, 0},
   { 0, 2,-1, 2, 0, 0},
   { 0, 0, 3, 0, 0, 0},
   {-2, 0, 0, 1, 1, 0},
   {-1, 0, 0,-1, 2,-3},
   {-1,-1, 0, 0, 0, 6}};

NAG represents a sparse matrix by the position and value of the non-zero elements. We can use Mathematica to extract the positions and non-zero values of the above matrix.

Extract the positions (rowcols) of the non-zero elements of the above matrix

rowcols = Position[ a, _?(#=!=0&), {2}, Heads->False ];
{therows,thecols} = Transpose[rowcols];

Extract the corresponding non-zero values (thevals)

thevals = Apply[ a[[#1,#2]]&, rowcols, {1} ];

Factorise the sparse matrix

Factorising a sparse matrix may lead to a less sparse matrix. One must provide some additional space for the factorised matrix and this is documented in the NAG manual.

workspace=Table[0,{4*Length[thevals]}];
f01brf[ Length[thevals],           (* $NZ *)
        Join[thevals,workspace],   (* $A *)
        Join[therows,workspace],   (* $IRN *)
        Join[thecols,workspace] ]; (* $ICN *)

Having factorised A, solve for some right hand side b,

b={1, 2, 3, 4, 5, 6};
f04axf[b]

(* % - Inverse[N[a]].N[b] //Chop *)
a=.;b=.; rowcols=.;thevals=.;therows=.;thecols=.; workspace=.

Solve Ax=b for A a very big sparse matrix

Here is a random 1000x1000 sparse matrix

We shall make a random matrix with roughly 3 non-zero elements per row.

n=1000; nz=3*n;
rcdiags = Transpose[{Range[n],Range[n]}];
rowcols = Table[ {Random[Integer,{1,n}],
                  Random[Integer,{1,n}]}, {nz-n} ];
rowcols = Join[ rcdiags, rowcols ]//Union; rcdiags=.;
{therows,thecols} = Transpose[rowcols];
thevals = Table[ Random[], {Length[rowcols]} ];

rowcols//Short

thevals//Short

Factorise the sparse matrix

The following takes about 110secs of Mathematica time, plus 30secs of NAG time.

workspace=Table[0,{4*Length[thevals]}];
f01brf[ Length[thevals],           (* $NZ *)
        Join[thevals,workspace],   (* $A *)
        Join[therows,workspace],   (* $IRN *)
        Join[thecols,workspace] ]; (* $ICN *)

Having factorised A, solve for some right hand side b,

The following takes about 3sec of Mathematica time.

b=Table[Random[],{n}];
x=f04axf[b]; x//Short

Check the solution

The following should return the zero vector if x is the solution. (This check takes about 110secs of Mathematica time to compute.)

Table[ Module[{p=Position[therows,i]//Flatten},
          thevals[[p]] . x[[thecols[[p]]]] - b[[i]] ],
       {i,1,n} ]//Short

n=.; nz=.; b=.; x=.
rowcols=.; thevals=.; therows=.; thecols=.; workspace=.

Similar IMSL Examples

find["dqdag " || "d01akf"]

find["dqdagi" || "d01amf"]

find["dqdags" || "d01ajf"]

find["dai " || "s17agf"]

find["dcscon" || "e01baf"]

find["dcsval" || "e02bbf"]

find["dneqnf" || "c05nbf"]

find["dzreal" || "c05adf"]

find["divprk" || "d02baf"]

find["???" || "d02ebf"]

find["dbvpfd" || "d02gaf"]

find["dlfcrg" || "f01btf" && Not["f04ayf"]]

find["dlfsrg" || "f04ayf"]

find["???" || "f04eaf"]

find["dlqrrr" || "f01qaf"]

find["devcsf" || "f02abf"]

find["devcrg" || "f02agf"]

find["???" || "f01mcf"]

find["???" || "f04mcf"]

find["???" || "f01brf"]

find["???" || "f04axf"]