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
------------------------------