matlab c混合编程_辛普森法则计算定积分

matlab C混合编程
Given function f(x) =sin(x)/x , the integral range [a, b], please work out its numericalintegral in range [a, b]. You are required to calculate the numerical integral by Simpsons rule.
使用辛普森法则计算函数f(x) =sin(x)/x 在[a,b]的定积分。

//myQuad.c
#include <mex.h>
#include <math.h>       //Introduce "math.h" header file
//#define N 10000000
int N = 0;
double dCalcuSumI(double x1, double x2);    
double dCalcuIntegral(double a, double b);

void mexFunction(int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[])      
{
    //Parameter name explanation:
    //rsh = right hand side, input parameters
    //lsh = left hand side, output parameters
    //prsh = pointer right hand side
    //nrsh = number right hand side
    //plsh = pointer left hand side
    //nlsh = number left hand side

    if(nrhs < 2)
        mexErrMsgTxt("Input arguments' num is wrong!n");       //Check the input arguments' number

    double a,b;

    a = mxGetScalar(prhs[0]);
    b = mxGetScalar(prhs[1]);
    //prhs[] is input parameter array, but the element isn't int/double... type, 
    //the element is a complex structure, if the input parameter is double type,
    //"mex.h" provide a function named "mxGetScalar(prhs[i])" to get the double value from this complex structure(mxArray)  
    N = (int)mxGetScalar(prhs[2]);
    plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL);
    //plhs[] is a empty array, so we should initialize the first element's type
    //output is also a double type data, so create a 1*1 matrix 
    //by using "mex.h"'s function "mxCreateDoubleMatrix(line num, row num, type)"
    double *dSum;
    dSum = mxGetPr(plhs[0]); 
    //"mxGetPr(plhs[i])" can get the first (double *) type pointer from plhs[]

    *dSum = dCalcuIntegral(a,b);
}

double dCalcuSumI(double x1, double x2)     //Calculate section sum
{
    double sum;
    sum = (x2 - x1) / 6 * (sin(x1)/x1 + sin(x2)/x2 + 4 * sin((x1+x2)/2) / ((x1+x2)/2));
    return sum;
}

double dCalcuIntegral(double a, double b)       //Calculate the Integral value via Simpson's rule
{
    double d = (b-a)/N;
    double dSum = 0;
    int i;
    for(i = 0; i < N; i++)
    {
        dSum += dCalcuSumI(a+i*d, a+(i+1)*d);
    }
    return dSum;
}

//test.m
clear;
interval = input('input the interval：')
interval = interval + 0.00001;
k = 8/interval;
for i = 1:k+1
    x(i) = -2 + i*interval;
    y(i) = myQuad(-2,-2 + i*interval);
end
plot(x,y,'-')