InterCall version 2 Examples Terry Robb 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 <>, D02EBF, D05ABF} InterCall::linked: Using remote driver version 2.0 on host VX23. 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) (* < Min[ROWS[$A], ROWS[$B]], (* DATA=I *) $A -> In, (* DATA=R[$NDIM] *) $B -> In, (* DATA=R[$NDIM] *) $MINPTS -> 1 :> d01fcf`minpts, (* DATA=I *) $MAXPTS -> 100, (* DATA=I *) $FUNCTN -> In, (* DATA=RF[I, R[$NDIM]] *) $EPS -> 0.00001, (* DATA=R *) $ACC :> d01fcf`acc, (* DATA=R *) $LENWRK -> ($NDIM + 2)*(1 + $MAXPTS/(2^$NDIM + 2*$NDIM^2 + 2*$NDIM + 1)), (* DATA=I *) $WRKSTR :> Null, (* DATA=R[$LENWRK] *) $FINVAL :> Out, (* DATA=R *) $IFAIL -> -1 (* DATA=I *) ] (* CODE="LIBRARY" *) d01fcf[$A_, $B_, $FUNCTN_] -> $FINVAL 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 d01fcf[$A_, $B_, $FUNCTN_] -> $FINVAL 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) *) 0.7641752396967671 Now we can compare the Mathematica and NAG answers, with the exact answer: {ans1,ans2} - N[2/3 + (Log[10]-CosIntegral[10]+EulerGamma)/30] -10 -14 {-6.15552 10 , 9.92539 10 } 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) *) 0.00698908 find["quadra" && "finit" && "one" && "oscil" && "(NAG)"] D01AKF (NAG) One-dimensional quadrature, adaptive integration over a finite interval, method suitable for oscillating functions. ans2=d01akf[ Function[x, x Cos[30x] Sin[x]], 0, 2Pi ] (* ShowTime::timing: 0.333333 Second (0.92579) *) 0.00698908265537221 {ans1,ans2} - N[2Pi/899] -14 -17 {-2.81112 10 , 3.46945 10 } 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) *) 3.14159 find[ "quadrature" && "infinite" && "one" && Not["singular"] && "(NAG)" ] D01AMF (NAG) One-dimensional quadrature, adaptive integration over an infinite or semi-infinite interval. ans2=d01amf[ Function[x, 1/((x+1)Sqrt[x])], 0, +1, $EPSREL->10^-14 ] (* ShowTime::timing: 0.383333 Second (0.9924) *) 3.14159265358979 {ans1,ans2} - N[Pi] -9 -15 {-6.70358 10 , -3.10862 10 } 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) *) -2.54326 find["one" && "quadr" && "finit" && "bad" && "(NAG)"] D01AJF (NAG) One-dimensional quadrature, adaptive integration over a finite interval, strategy due to Piessens & deDoncker, allowing for badly-behaved integrands. ans2=d01ajf[ Function[x,x Sin[30x]/Sqrt[1-(x/(2Pi))^2]], 0, 2Pi, $EPSREL->10^-9 ] (* ShowTime::timing: 0.533333 Second (1.34366) *) -2.54325961889401 {ans1,ans2} - N[2*Pi^3*BesselJ[1, 60*Pi]] -11 -13 {1.40146 10 , -4.23661 10 } 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) - l º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)"] D05AAF (NAG) Linear non-singular Fredholm integral equation, 2nd kind, split kernel. D05ABF (NAG) Linear non-singular Fredholm integral equation, 2nd kind, smooth kernel. 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 f[x_] := 1. - 1.499999999999396*x + 0.4999999999592929*x^2 - 0.1666666658533438*x^3 + 0.04166665915372853*x^4 - 0.00833329486912761*x^5 + 0.001388769492590569*x^6 - 0.000198178094865398*x^7 + 0.00002450662868795917*x^8 - 2.52334093602258*10^-6*x^9 + 1.680948116700167*10^-7*x^10 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 -14 1.28786 10 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)"] F01BRF (NAG) LU-factorisation, real sparse matrix. F01BSF (NAG) LU-factorisation, real sparse matrix with known sparsity pattern. F04AXF (NAG) Simultaneous linear equations, (coefficient matrix already factorised), real sparse matrix, approximate solution. F04QAF (NAG) Least-squares, m real equations in n unknowns, sparse linear least-squares problem. 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} ] {5, 2, -1, 2, 3, -2, 1, 1, -1, -1, 2, -3, -1, -1, 6} 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] {0.2, 1.6, 1., -0.1, 4.5, 1.3} And the above is the solution x. As a check we could multiply it by A and compare with b, a.% - b //Chop {0, 0, 0, 0, 0, 0} 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)" ] D02BAF (NAG) System of ordinary differential equations, initial value problem, Runge-Kutta-Merson method, (simple driver), over a range. D02BBF (NAG) System of ordinary differential equations, initial value problem, Runge-Kutta-Merson method, (simple driver), over a range with intermediate output. D02BGF (NAG) System of ordinary differential equations, initial value problem, Runge-Kutta-Merson method, (simple driver), until a component of the solution attains a given value. D02BHF (NAG) System of ordinary differential equations, initial value problem, Runge-Kutta-Merson method, (simple driver), until a function of the solution is zero. 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.