MATLAB 与 C 语言的混合编程实战之辛普森积分法、自适应辛普森积分

//Author：glm
#include<cstdio>
#include<cmath>
#include<iostream>
#include "mex.h"
#define ll long long int
#define rg register ll

inline double f(double x)
{
    if(x==0)return 1;
    return sin(x)/x;
}
inline double calculate(double a,double b)//int(f,a,b)=(b-a)/6*(f(a)+4*f((a+b)/2+f(b))
{
    double sum=0;
    for(rg i=0;i<(ll)((b-a)/0.025);i++)
    {
        double a1=a+0.0250*i,b1=a1+0.0250,mid=(a1+b1)/2.0;
        sum+=(mid-a1)/6*(f(a1)+4*f((a1+mid)/2)+f(mid))+(b1-mid)/6*(f(mid)+4*f((b1+mid)/2)+f(b1));
        
    }
    return sum;
   /* int parts=(int)((b-a)/0.01);
    double ans=0;
    for (int i=0;i<parts;i++)
        {
            double ai=a+i*0.01;
            double bi=ai+0.01;
            double mi=(ai+bi)/2;
            ans+=(mi-ai)/6*(f(ai)+4*f((ai+mi)/2)+f(mi))+(bi-mi)/6*(f(mi)+4*f((bi+mi)/2)+f(bi));
            //ans+=(bi-ai)/6*(f(ai)+4*f((ai+bi)/2)+f(bi));
        }
    return ans;*/
}
void mexFunction (int nlhs,mxArray *plhs[],int nrhs,const mxArray * prhs[])
{
    double *a;
    double b,c;
    plhs[0]=mxCreateDoubleMatrix(1,1,mxREAL);
    a=mxGetPr(plhs[0]);//
    b=*(mxGetPr(prhs[0]));
    c=*(mxGetPr(prhs[1]));
    *a=calculate(b,c);
}

documentclass[12pt]{article}
usepackage{amsmath}
usepackage{mathtools}
usepackage{graphics,epsfig}
usepackage{xspace}
usepackage{color}
usepackage{amssymb}
usepackage{subfigure}
usepackage{multirow}
usepackage{listings}
usepackage{graphicx}
usepackage{url}
%%uncomment following line if you are going to use ``Chinese''
%usepackage{ctex} 

definecolor{ballblue}{rgb}{0.13, 0.67, 0.8}
definecolor{cornflowerblue}{rgb}{0.39,0.58,0.93}
definecolor{babyblueeyes}{rgb}{0.63, 0.79, 0.95}

% preset-listing options
lstset{
  backgroundcolor=color{white},   
  basicstyle=footnotesize,    
  language=matlab,
  breakatwhitespace=false,         
  breaklines=true,                 % sets automatic line breaking
  captionpos=b,                    % sets the caption-position to bottom
  commentstyle=color{ballblue},    % comment style
  extendedchars=true,              
  frame=single,                    % adds a frame around the code
  keepspaces=true,                 
  keywordstyle=color{blue},       % keyword style
  numbers=left,                    
  numbersep=5pt,                   
  numberstyle=tinycolor{blue}, 
  rulecolor=color{babyblueeyes},
  stepnumber=1,              
  stringstyle=color{black},     % string literal style
  tabsize=4,                       % sets default tabsize to 4 spaces
  title=lstname                   
}

usepackage{geometry}
geometry{
 a4paper,
 total={210mm,297mm},
 left=20mm,
 right=20mm,
 top=20mm,
 bottom=20mm,
 }

marginparwidth = 10pt


begin{document}
title{Weekly Assignment:Simpson’s rule of calculation}
author{ \ Stdudent  Guo LiMin \ Id: 22920182204174}
maketitle


section{Problem Description}

Given function  $f(x)=cfrac{sinx}{x}$ , the integral range [a, b], please work out its numerical integral in range [a, b]. 
You are required to calculate the numerical integral by Simpson’s rule.
\\{color{red}Requirements:}\\    1. Your function myQuad MUST be implemented in C and called by Matlab\
2. The command interface in Matlab looks like: v = myQuad(a, b);\
3. You are NOT allowed to use “recursive procedure” when you implement Simpson’s
rule\
4. Given $F(x) =int^{6}_{-2}cfrac{sinx}{x}dx$, please plot out the curve of F(x) in the range [-2 6]\
5. According to Simpson’s rule, the number of intervals will impact on the precision.
You are required to try different numbers of intervals and see how it impacts your
result. Please present an analysis about this in your report.

section{Code of Solution}
The followings are codes implemented by C++,no fancy algorithm included,working out this problem only by using simple simulations\
%
%if one wants to paste part of the code into the report
%one can put in following manner
begin{lstlisting}[language=c++, linewidth=linewidth,caption={main working function(myQuad.cpp)}]
#include<cstdio>
#include<cmath>
#include<iostream>
#include "mex.h"
#define ll long long int
#define rg register ll
inline double f(double x)
{
    return sin(x)/x;
}
inline double calculate(double a,double b)//int(f,a,b)=(b-a)/6*(f(a)+4*f((a+b)/2+f(b))
{
    double sum=0;
    for(rg i=0;i<(ll)((b-a)/0.001);i++)
    {
        double a1=a+0.001*i,b1=a1+0.001,mid=(a1+b1)/2.0;
        sum+=(mid-a1)/6*(f(a1)+4*f((a1+mid)/2)+f(mid))+(b1-mid)/6*(f(mid)+4*f((b1+mid)/2)+f(b1));
        
    }
    return sum;
   /* int parts=(int)((b-a)/0.01);
    double ans=0;
    for (int i=0;i<parts;i++)
        {
            double ai=a+i*0.01;
            double bi=ai+0.01;
            double mi=(ai+bi)/2;
            ans+=(mi-ai)/6*(f(ai)+4*f((ai+mi)/2)+f(mi))+(bi-mi)/6*(f(mi)+4*f((bi+mi)/2)+f(bi));
            //ans+=(bi-ai)/6*(f(ai)+4*f((ai+bi)/2)+f(bi));
        }
    return ans;*/
}
void mexFunction (int nlhs,mxArray *plhs[],int nrhs,const mxArray * prhs[])
{
    double *a;
    double b,c;
    plhs[0]=mxCreateDoubleMatrix(1,1,mxREAL);
    a=mxGetPr(plhs[0]);//
    b=*(mxGetPr(prhs[0]));
    c=*(mxGetPr(prhs[1]));
    *a=calculate(b,c);
}

end{lstlisting}

begin{lstlisting}[language=c++, linewidth=linewidth,caption={Draw the contrast curve and main interface code(test.m)}]
mex myQuad.cpp
syms t
a=-2:0.27:6;
b=-2:0.27:6;
c=-2:0.27:6;
cnt=0;
for x=-2:0.27:6
y=int(sin(t)/t,-2,x);
cnt=cnt+1;
b(cnt)=y;
c(cnt)=myQuad(-2,x);
b(cnt),c(cnt);
end
plot(a,b,'r',a,c,'b')
    
end{lstlisting}
section{Experiment Theory and Results}
Given function f(x), and its range of x [a,b], it is possible to estimate its integral in range [a,b]
by quandratic interpolation based on “Simpson’s rule”. The “Simpson’s rule” is given as
follows.\
begin{equation}
begin{split}$int^{b}_{a}cfrac{sinx}{x}dx &=cfrac{(b-a)}{6}*[f(a)+4*f(cfrac{a+b}{2})+f(b)]\
            &=cfrac{(m-6)}{6}*[f(a)+4*f(cfrac{a+m}{2})+f(m)]+cfrac{(b-m)}{6}*[f(m)+4*f(cfrac{b+m}{2})+f(b)]nonumber$   
end{split}
end{equation}
\
What’s more,The more number of segments used in the range [a, b], the better is the approximation.
\\
To start with,we should learn how to compile c/c++ files with matlab and how to write fantastic and unified-formatted papers,
this aspect of learning has been updated to my personal blog  {color{red}url{newcoder-glm.blog.csdn.net}},so for more information,please go up to visit my blog.
\\
After mastering the mixtrue of C/C++& matlab programming technique(It's just an interface that calls an external compiler in essence),as easy as pie we can write the codes above and 
should curve figures(check out the attachment for more details) which shows graphs of the result and the exact result obtained by Simpson integral method when the interval of integral interval is divided differently(when the value is 0.05 0.25)

begin{table}[h]
caption{Circumstances when we apply different intervals}
begin{center}
	begin{tabular}{|c|c|c|c|}
	hline
	Interval & 0.01 & 0.10 & 0.25\ hline
	Results & textbf{2.551496047169967}& textbf{2.551496048068467} & textbf{2.551496047173387  } \ hline
	end{tabular}
end{center}
end{table}


subsection{further ideas}
We can see from the Simpson integral (the simple three-point formula) that there is some error, and it's easy to see that the smaller the interval, the more accurate the integral,
however,if we shorten the interval, of course the results' accuracy come up, but at the cost of higher time complexity.\
So here we introduce "Adaptive Simpson’s rule" which is the integral can automatically control the size and length of the cutting interval, so that the accuracy can meet the requirements.
In particular, if the midpoint of [a,b] is c, the result will be directly returned if and only then {color{red}$[S(a,c)+S(c,b)-S(a,b)]<15Eps$}  Otherwise, the recursive call will divide the interval again. The accuracy of Eps should be reduced by half when recursion is called.
begin{figure}[h]
	centering
	{includegraphics[width=0.7linewidth]{025.pdf}}
    hspace{0.15in}
end{figure}

begin{figure}[h]
	centering
	{includegraphics[width=0.7linewidth]{005.pdf}}
    hspace{0.15in}
end{figure}

section{What I learned from the Project}
1.learning how to programing cpp files with powerful matlab,excuted with high efficiency and conciseness}
\\2.the ability of exploring new territory that never have been leart and mastering it(hand in report in time)}
\\3.learning how to write fantastic and unified-formatted papers with modernized Latex}

section{Two attached pictures of point three}
end{document}