Mathematica 9 is now available
Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2010

[Date Index] [Thread Index] [Author Index]

Search the Archive

Mathlink problem in C

  • To: mathgroup at smc.vnet.net
  • Subject: [mg112214] Mathlink problem in C
  • From: rsgill <rsgill at phas.ubc.ca>
  • Date: Sun, 5 Sep 2010 05:26:45 -0400 (EDT)

Hi,

I am trying to learn mathlink and how it can be used with C
programming. To that end ,
I am attempting to write a gaussian quadrature scheme in C to do
integrals numerically. Also, I am using Mathematica's
gaussianquadratureweights function which gives you the weights and
abcissas for the integral approximation.

I am well aware of the fact that routines in C are available to do
this, but I am using
quadrature as a test case to learn mathlink.

In order to use the function "gaussianquadratureweights", one has to
first load the package "NumericalDifferentialEquationAnalysis".

[Problem]: The problem is that I can't get mathlink to give me back
the
right list of weights and abcissas for the above function.

[Code]: In my C code I have defined a function called "gaussquad"
which
takes in the number of intervals n, xmin value a, and xmax value b for
integration of some function , say sin(x) over the interval [a,b].
Also,
my gaussquad should return two vectors xi[n] and ci[n] to the main
program and these are the abcissas and weights, respectively.

Then, one can simply do a sum, essentially a dot product of the two
vectors, in the main program to obtain an approximation to the
integral.

I am attaching my C code and also the mathematica file to do gaussian
quadrature.

--------------------C code ----------------------------------------
/*Gaussian Quadrature scheme for integration*/

#include <stdio.h>
#include <math.h>
#include <mathlink.h>

void gaussquad(MLINK lp, int, double, double, double *, double *);
double func(double);
int main()
{
        int n, i;
        double a, b, sum;
/*--------------Mathematica Intialization-------------------*/
        int argc = 4;
        char *argv[5] = {"-linkname", "math -mathlink", NULL, NULL,
NULL};
        MLINK lp;
        int pkt;
        MLEnvironment env;

        env = MLInitialize(NULL);
        if (env == NULL) return 1;
        lp = MLOpen(argc, argv);
        if (lp == NULL) return 1;
/*--------Quadrature Scheme---------------*/
        scanf("%d %lf %lf", &n, &a, &b);   //Input order n, and
integration limits
        double ci[n], xi[n];    //Defines two vectors: weights and
abcissas

        gaussquad(lp, n, a, b, (double *)ci, (double *)xi);   //Send
n, a, b to mathematica function,
                                                                                //
to get ci and xi vectors
        for (i = 0, sum = 0; i < n; i++) {        //Perform
approximation to integral
                sum += ci[i]*func(xi[i]);
        }
        printf("%lf\n", sum);
/*--------Closing Mathematica Link-------------------*/
        MLPutFunction(lp, "Exit", 0);
        MLEndPacket(lp);
        MLClose(lp);
        MLDeinitialize(env);
        return 0;
/
*--------------------------------------------------------------*/
}

double func(double x)    //Evaluates any given function at abcissas
from
                                   //gaussquad
{
        double funcval;

        funcval = sqrt(1 + x*x);
        return funcval;
}

void gaussquad(MLINK lp, int n, double a, double b, double *x, double
*c)
{
        long length, argcount;
        const char *funcname;
        int intarg;


        MLPutFunction(lp, "EvaluatePacket", 1);
                MLPutFunction(lp, "Get", 1);
                        MLPutString(lp,
"\"NumericalDifferentialEquationAnalysis`\"");
        MLEndPacket(lp);
        MLNewPacket(lp);
        MLPutFunction(lp, "EvaluatePacket", 1);
                MLPutFunction(lp, "Set", 2);
                        MLPutSymbol(lp, "coord");
                        MLPutFunction(lp, "GaussianQuadratureWeights",
3);
                                MLPutInteger(lp, n);
                                MLPutDouble(lp, a);
                                MLPutDouble(lp, b);
        MLEndPacket(lp);
        printf("Func Type: %d, %s, %d\n", MLGetFunction(lp, &funcname,
&intarg), funcname, intarg);
        printf("Return: %d, %ld\n", MLCheckFunction(lp,
"ReturnPacket",
&argcount), argcount);
        MLNewPacket(lp);
        MLPutFunction(lp, "EvaluatePacket", 1);
                MLPutFunction(lp, "Part", 3);
                        MLPutSymbol(lp, "coord");
                        MLPutString(lp, "All");
                        MLPutInteger(lp, 1);
        MLEndPacket(lp);
        printf("Return: %d, %ld\n", MLCheckFunction(lp,
"ReturnPacket",
&argcount), argcount);
        while (MLNextPacket(lp) != RETURNPKT) MLNewPacket(lp);
        MLGetRealList(lp, &x, &length);
        printf("Error: %d\n", MLError(lp));
        printf("%s\n", MLErrorMessage(lp));
        MLNewPacket(lp);
        MLPutFunction(lp, "EvaluatePacket", 1);
                MLPutFunction(lp, "Part", 3);
                        MLPutSymbol(lp, "coord");
                        MLPutString(lp, "All");
                        MLPutInteger(lp, 2);
        MLEndPacket(lp);
        printf("Return: %d, %ld\n", MLCheckFunction(lp,
"ReturnPacket",
&argcount), argcount);
        while (MLNextPacket(lp) != RETURNPKT) MLNewPacket(lp);
        MLGetRealList(lp, &c, &length);
}

-------------------Mathematica
Code------------------------------------------------
Remove["Global`*"]
<<NumericalDifferentialEquationAnalysis`
F[x_]:=Sqrt[1+x^2]
Int[n_,a_,b_]:=Module[{i,coord,x,w},coord=GaussianQuadratureWeights[n,
0,2];Table[Subscript[x, i]=coord[[i]][[1]],{i,n}];Table[Subscript[w,
i]=coord[[i]][[2]],{i,n}];Sum[Subscript[w, i]*F[Subscript[x, i]],
{i,n}]]
Int[6,0,2]
2.95789
NIntegrate[F[y],{y,0,2}]
2.95789
------------------------------


  • Prev by Date: Re: FindRoots?
  • Next by Date: Re: FindRoots?
  • Previous by thread: Re: NDSolve -- indexing of dependent variable that is arbitrary list
  • Next by thread: Finite Groups...infinite disappoinment