• 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

In order to use the function "gaussianquadratureweights", one has to

[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

--------------------C code ----------------------------------------

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

double func(double);
int main()
{
int n, i;
double a, b, sum;
/*--------------Mathematica Intialization-------------------*/
int argc = 4;
NULL};
int pkt;
MLEnvironment env;

env = MLInitialize(NULL);
if (env == NULL) return 1;
lp = MLOpen(argc, argv);
if (lp == NULL) return 1;
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);
MLPutFunction(lp, "Exit", 0);
MLEndPacket(lp);
MLClose(lp);
MLDeinitialize(env);
return 0;
/
*--------------------------------------------------------------*/
}

double func(double x)    //Evaluates any given function at abcissas
from
{
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");
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]
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