InterCall version 2

Examples

Terry Robb
<tdr@vaxc.cc.monash.edu.au>

Abstract

This Notebook contains some examples which require the InterCall compiler available in version 2 of InterCall. It includes examples in the following areas:
   – Quadrature (one-dimensional and multi-dimensional)
   – Integral equations
   – Sparse matrices
   – Ordinary differential equations (including stiff systems)

In most of the examples, InterCall easily wins in real-time against equivalent Mathematica routines. When it loses it is usually because the compilation of a passed function is slow. This is because the InterCall compiler is currently written as a Mathematica routine. Eventually, however the compiler will be written in C to speed up compilation.

These examples require the NAG subroutine library. There are several other example Notebooks, available via MathSource, that use public domain libraries and require the InterCall compiler. A Notebook solving the Traveling Salesman Problem can be obtained by emailing the line
        find 0203-432
to mathsource@wri.com

Initialisation

Load in the InterCall package and NAG database

<<InterCall.m; <<InterData.m;

Import in several routines from a site that has NAG

The NAG library we are going to access is not on the same computer as the Mathematica kernel. Here is the name of a computer that has NAG,

myNagSite= "vaxc.cc.monash.edu.au"; (* EDIT THIS LINE *)

From that computer we will eventually import the following NAG routines,

routines = {D01AKF,D01AMF,D01AJF,D01FCF,D02BBF,D02EBF,D05ABF};

To connect to the computer with NAG we use the InterCall command, and to import in some routines we use Import.

imported = (InterCall[myNagSite]; Import[routines]);

If for some reason we couldn't import all the routines, then we shall quit now.

If[ imported===Null, Quit[] ];
myNagSite=.; routines=.; imported=.;

Load in a Timing package (recording cpu and elapsed time)

(* <<MyShowTime.m; *)

Quadrature

Multi-dimensional quadrature

Integrate x2+y2+sin(10xy)/3 over the unit square using Mathematica then NAG. The exact answer is 2/3 + (ln(10)-Ci(10) + Γ)/30.

Here is a plot of the function over the region to be integrated.

Plot3D[ x*x+y*y+Sin[10 x y]/3, {x,0,1}, {y,0,1} ];

We could numerically integrate it with Mathematica's NIntegrate command as follows:

ans1 = NIntegrate[x*x+y*y+Sin[10 x y]/3, {x,0,1},{y,0,1}]
(* ShowTime::timing: 7.38333 Second (8.20319) *)

But suppose we wanted to use a NAG routine as a cross-check. First we need to find an appropriate NAG routine that can do adaptive multi-dimensional quadrature. InterCall's find command can be used to search for such a routine:

find[ "multi" && "quadratu" && "adapt" && "(NAG)" ]

There are thus two possible routines, and we shall choose d01fcf. In order to call that routine correctly we must read the NAG manual and also use InterCall's GetDefault command to see what defaults are set and also to see what the arguments mean.

GetDefault[d01fcf]

The above shows quite a lot of detail about the default value and datatype of each argument to the routine d01fcf. The final output line is however the main piece of information to focus on, and this line alone could be obtained by simply typing

?d01fcf

By reading the NAG manual, we see that A and B are the limits of the integration and that FUNCTN is a function defining the integrand. The manual also says that FINVAL is the integral. An appropriate call to d01fcf is thus:

ans2 = d01fcf[ {0,0}, {1,1},
               Function[{n,xn},Module[{x,y},{x,y}=xn;
                               x*x+y*y+Sin[10 x y]/3]],
               $EPS->10^-9, $MAXPTS->30000 ]
(* ShowTime::timing: 0.733333 Second (3.44138) *)

Now we can compare the Mathematica and NAG answers, with the exact answer:

{ans1,ans2} - N[2/3 +
              (Log[10]-CosIntegral[10]+EulerGamma)/30]

Overall we see that calling the NAG routine using InterCall is quicker in real-time than calling Mathematica's NIntegrate routine.

ans1=.; ans2=.;

Quadrature of oscillating function over finite interval

Integrate x cos(30x) sin(20x) from 0 to 2π using Mathematica then NAG. The exact answer is 2π / 899.

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

ans1=NIntegrate[ x Cos[30x] Sin[x], {x, 0, 2Pi} ]
(* ShowTime::timing: 5.21667 Second (5.40272) *)

find["quadra" && "finit" && "one" && "oscil" && "(NAG)"]

ans2=d01akf[ Function[x, x Cos[30x] Sin[x]], 0, 2Pi ]
(* ShowTime::timing: 0.333333 Second (0.92579) *)

{ans1,ans2} - N[2Pi/899]

ans1=.; ans2=.;

Quadrature over semi-infinite interval

Integrate 1/((x+1)√x) from 0 to Infinity using Mathematica then NAG. The exact answer is π.

Plot[ 1/((x+1)Sqrt[x]), {x,0.01,10} ];

ans1=NIntegrate[ 1/((x+1)Sqrt[x]), {x, 0, Infinity} ]
(* ShowTime::timing: 1.21667 Second (1.33959) *)

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

ans2=d01amf[ Function[x, 1/((x+1)Sqrt[x])], 0, +1,
             $EPSREL->10^-14 ]
(* ShowTime::timing: 0.383333 Second (0.9924) *)

{ans1,ans2} - N[Pi]

ans1=.; ans2=.;

Quadrature of badly-behaved integrand over finite interval

Integrate x sin(30x)/√(1-(x/2π)2) from 0 to 2π using Mathematica then NAG. The exact answer is 2 π3 J1(60π).

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

ans1=NIntegrate[ x Sin[30x]/Sqrt[1-(x/(2Pi))^2],
                {x, 0, 2Pi-10^(-15)} ]
(* ShowTime::timing: 4.06667 Second (4.64816) *)

find["one" && "quadr" && "finit" && "bad" && "(NAG)"]

ans2=d01ajf[ Function[x,x Sin[30x]/Sqrt[1-(x/(2Pi))^2]],
             0, 2Pi, $EPSREL->10^-9 ]
(* ShowTime::timing: 0.533333 Second (1.34366) *)

{ans1,ans2} - N[2*Pi^3*BesselJ[1, 60*Pi]]

ans1=.; ans2=.;

Integral Equations

Fredholm integral equation, 2nd kind, smooth kernel

Solve the Fredholm integral equation using NAG. The equation to be solved for the function f(x) is f(x) - λ ∫ab k(x,t) f(t) dt = g(x), with k(x,t) and g(x) etc given as follows:

( k=Function[{x,t}, x Exp[t]];
  g=Function[x, Exp[-x]];
  lambda = -1; a=0; b=1; )

First we try find a NAG routine that can solve this type of integral equation.

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

We choose the routine d05abf. The following returns the coefficients for a eleventh order Chebyshev expansion of the solution to f(x) + ∫01 x et f(t) dt = e-x . The solution f(x) is neither an odd nor an even function, so we set $ODOREV to False.

c = d05abf[k,g,lambda,a,b,$ODOREV->False,$N->11];
(* ShowTime::timing: 0.666667 Second (1.07421) *)

Having found the Chebyshev coefficents c we can use them to construct the solution f(x) as follows:

f[x_] := Evaluate @ Expand[ First[c]/2 +
        Rest[c].Table[ChebyshevT[i,2x-1],
                      {i,Length[Rest[c]]}] ];

Here is the expansion for f(x).

?f

The analytic solution is e-x - x/2 and we can use it to check the accuracy of the numerical solution. A plot of the difference shows that the numerical solution is accurate to 1 part in 1014 over the interval [0,1].

Plot[ ( (Exp[-x]-x/2) - f[x] )* 10^14, {x,0,1} ];

Table[ (Exp[-x]-x/2) - f[x], {x,0,1,0.05} ]//Abs//Max

lambda=.; a=.; b=.; k=.; g=.; c=.; ClearAll[f]

Sparse matrices

Large, but sparse, matrices often occur in many problems including the solution of elliptic partial differential equations to name just one. Sparse matrices are an area in which Mathematica is currently deficient, although they can be handled relatively easily by some standard numerical libraries.

Solve Ax=b for A a sparse matrix

Suppose we have reduced a problem, by using Mathematica, down to solving real non-symmetric sparse matrix. By using InterCall we can find all relevant library routines that may be appropriate.

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

Two of the routines above are appropriate for this example – f01brf and f04axf. We can inspect their default setting, but for now we just assume we know how to use them and continue on.

Suppose the sparse system involves the following matrix (which has lots of zero's).

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 such a sparse matrix by the position and value of the non-zero elements. We can use Mathematica to extract those positions and non-zero values by,

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

The above finds the row and column number of any non-zero element in the matrix, and with this information we can go back and get the non-zero elements (thevals) of the matrix.

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

According to the NAG manual the next step is to factorise the matrix using f01brf. Factorising a sparse matrix may lead to a less sparse matrix, and one must provide some additional space for the factorised matrix. This is documented in the NAG manual. Anyway the incantation that will factorise any sparse matrix A is this:

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, it is now possible to solve the system for some specific right hand side b. The NAG manual suggests f04axf be used, thus

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

And the above is the solution x. As a check we could multiply it by A and compare with b,

a.% - b //Chop

Since the above is all zeroes we confirm that the answer is correct.

Bigger matrices

Although the system was just a 6x6 sparse matrix, the method can be applied to vastly bigger problems. Even a 1000x1000 sparse matrix is still fairly trivial and can be easily handled.

Ordinary Differential Equations

System of ODEs, initial value problem

Solve the Lorentz attractor equations using Mathematica then NAG. NDSolve beats d02bbf by a few seconds here only because some passed subroutines are slow to compile.

Mathematica method (NDSolve)

doit1[a_] := ( ans1 = NDSolve[
 { y1'[x] == -3 (y1[x] - y2[x]),            y1[0] == 0,
   y2'[x] == -y1[x] y3[x]- y2[x] + a y1[x], y2[0] == 1,
   y3'[x] ==  y1[x] y2[x]- y3[x],           y3[0] == 0 },
 { y1, y2, y3 }, {x,0,20}, MaxSteps->3000 ]; )

doit1[26.5];
(* ShowTime::timing: 3.46667 Second (3.69139) *)

ParametricPlot3D[ Evaluate[{y1[x],y2[x],y3[x]} /. ans1],
                  {x,0,20}, PlotPoints->1000 ];
(* ShowTime::timing: 18. Second (35.08594) *)

ans1=.; ClearAll[doit1];

NAG method (d02bbf)

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

Most NAG ODE routines have quite an awkward calling syntax from Mathematica's point of view. In the following we have to fight quite hard to call d02bbf. But it is nowhere near as painful as trying to call the routine from a fortran program. Although calling d02bbf using InterCall is not as easy as calling Mathematica's NDSolve it does illustrate the point that it is possible to access quite difficult external routines.

doit2[a_] := Block[{x0=0, y0={0,1,0}, xend=20,
                     mysub, myout, dxout, soln},
        mysub = Function[{x,y,f},
           Module[{y1,y2,y3}, {y1,y2,y3}=y;
                f = {-3 (y1 - y2),
                     -y1 y3 - y2 + a y1,
                      y1 y2 - y3} ]];
        myout = Function[{x,y},
           Module[{y1,y2,y3}, {y1,y2,y3}=y;
                AppendTo[soln, y1];
                AppendTo[soln, y2];
                AppendTo[soln, y3];
                x += dxout ]];
        soln={}; dxout=(xend-x0)/1000.;
        d02bbf[x0,xend,y0,mysub,myout];
        ans2 = Partition[soln,3]; ]

doit2[26.5];
(* ShowTime::timing: 2.9 Second (7.57823) *)

Show[ Graphics3D[ {Line[ans2]} ], Axes->Automatic ];
(* ShowTime::timing: 2.36667 Second (15.09379) *)

ans2=.; ClearAll[doit2];

Stiff system of ODEs, initial value problem

Solve a stiff system of ODEs using NAG.

( mysub = Function[{x,y,f}, Module[{f1,f2,f3,y1,y2,y3},
                {y1,y2,y3} = y;
                f1 = -0.04 y1 + 10000. y2 y3;
                f2 = -f1 - 30000000. y2 y2;
                f3 =       30000000. y2 y2;
                f = {f1,f2,f3}; ]];
  myder = Function[{x,y,pw}, Module[{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 }}; ]];
  myout = Function[{x,y},
                AppendTo[y1, y[[1]]];
                AppendTo[y2, y[[2]]];
                AppendTo[y3, y[[3]]];
                x = xnodes[[node]]; node++ ]; )

( xnodes = Table[ 10^(-6. + 0.1x), {x,0,71} ];
  xvals = Prepend[Drop[xnodes,-1], 0]; )

Most of the cpu-time used here is to compile mysub, myder and myout

( y1={}; y2={}; y3={}; node=1;
  x0 = First[xvals]; y0 = {1,0,0}; xend = Last[xvals];
  d02ebf[ x0, xend, y0, mysub, myder, myout,
          $IRELAB->2, $TOL->0.0001 ]; )
(* ShowTime::timing: 4.33333 Second (7.97996) *)

ListPlot[ Transpose[{xvals,y1}], PlotJoined->True ];
ListPlot[ Transpose[{xvals,y2}], PlotJoined->True ];
ListPlot[ Transpose[{xvals,y3}], PlotJoined->True ];
ListPlot[ Transpose[{Drop[xvals,1]//Log,Drop[y2,1]}],
                                 PlotJoined->True ];

( mysub=.; myder=.; myout=.; y1=.;y2=.;y3=.; node=.;
  xnodes=.; xvals=.; x0=.; y0=.; xend=.; )

References

InterCall

InterCall was developed by Terry Robb. For more information on InterCall send email to intercall_forward@wri.com. An order form is available via MathSource – email the line send 0202-587-0011 to mathsource@wri.com for details. Information can also be posted or faxed to Analytica below.

Analytica                     Fax: +61 (9) 386 5666
P.O. Box 522
Nedlands 6009
Perth WA
AUSTRALIA.