【工程数学】若干种解定积分的算法

// ConsoleAppDefInteSolu.cpp : 定义控制台应用程序的入口点。
//
/*
*函数功能：梯形公式与辛普生公式以及四阶Cotes公式直接求解定积分，这三种算法均是牛顿科次求积公式的1,2,4阶形式
*函数原形：
*double TrapezoidSolu(double a, double b),
*double SimpsonSolu(double a,double b);
*double Cotes4Solu(double a,double b);
*参数：double a, double b,积分区间的下限和上限
*返回值：积分结果
*时间复杂度：O(1)
*备注：三种简单算法的误差均较大，特别是梯形公式法
*日期：2014/11/22
*原创：是
*作者：EbowTang
*Email：tangyibiao520@163.com
*/

#include "stdafx.h"
#include "iostream"
using namespace std;

const double pi=3.1415926535897931;
double TrapezoidSolu(double a,double b);
double SimpsonSolu(double a,double b);
double Cotes4Solu(double a,double b);
int _tmain(int argc, _TCHAR* argv[])
{
double r0=0.0,r1=0.0,r2=0.0;//保存结果
double a=1.0;
double b=2.0;

cout<<"求函数“(exp(-x)-sin(pi/2))/x”在1到2上的积分"<<endl;

cout<<"/*********梯形公式的华丽分割线***********/"<<endl;
r0= TrapezoidSolu(a,b);
cout<<"梯形公式积分的结果为："<<r0<<endl;

cout<<"/*********辛普生公式的华丽分割线*********/"<<endl;
r1= SimpsonSolu(a,b);
cout<<"辛普生公式积分的结果为："<<r1<<endl;

cout<<"/********Cotes公式的华丽分割线***********/"<<endl;
r2=Cotes4Solu(a,b);
cout<<"科次4阶公式积分的结果为："<<r2<<endl;
system("pause");
return 0;
}

/*******被积函数***********/
double func(double x)
{
return (exp(-x)-sin(pi/2*x))/x;
}

/********辛普生计算函数*************/
double SimpsonSolu(double a, double b)
{
double h=(b-a)/2;
return (b-a)/6*(func(a)+func(b)+4*func(a+h));//辛普生公式直接求解（误差较大）
}

/********梯形计算函数*************/
double TrapezoidSolu(double a, double b)
{
return (b-a)/2*(func(a)+func(b));//梯形公式直接求解(误差大)
}

double Cotes4Solu(double a,double b)
{
double h=(b-a)/4;
return (b-a)/90*(7*func(a)+32*func(a+h)+12*func(a+2*h)+32*func(a+3*h)+7*func(b));//Cotes公式直接求解(误差较大)
}

// ConsoleAppCompositeTrapezoidSolu.cpp : 定义控制台应用程序的入口点。
//
/*
*函数功能：复化梯形公式以及复化辛普生公式求定积分
*函数原形：
*double CompositeTrapezoidSolu(double a, double b,int n),
*double CompositeSimpsonSolu(double a,double b,int n);
*参数：double a, double b,积分区间的下限和上限;int n,将要分割的段数
*返回值：积分结果
*时间复杂度：O(1)
*备注：误差相比于上面的三种算法将大大提高
*日期：2014/11/22
*原创：是
*作者：EbowTang
*Email：tangyibiao520@163.com
*/

#include "stdafx.h"
#include "iostream"
using namespace std;

const double pi=3.1415926535897931;

double CompositeTrapezoidSolu(double a,double b,int n);
double CompositeSimpsonSolu(double a,double b,int n);
double TrapezoidSolu(double a, double b);
double SimpsonSolu(double a, double b);

int _tmain(int argc, _TCHAR* argv[])
{
double r0=0.0,r1=0.0;//保存结果
double a=1.0;
double b=2.0;
int n=4;
char ans='n';
do
{
cout<<"求函数“(exp(-x)-sin(pi/2))/x”在1到2上的积分"<<endl;
cout<<"您想将积分区间分成多少段?"<<endl;
cin>>n;
cout<<"/*********复化梯形公式求积分的华丽分割线***********/"<<endl;
r0= CompositeTrapezoidSolu(a,b,n);
cout<<"复化梯形公式积分的结果为："<<r0<<endl;

cout<<"/*********复化辛普生公式求积分的华丽分割线*********/"<<endl;
r1= CompositeSimpsonSolu(a,b,n);
cout<<"复化辛普生公式积分的结果为："<<r1<<endl;
cout<<"是否继续？（y/n）:";
cin>>ans;
} while (ans=='y');
return 0;
}

/*******被积函数***********/
double func(double x)
{
return (exp(-x)-sin(pi/2*x))/x;
}

/********复化辛普生求积分计算函数*************/
double CompositeSimpsonSolu(double a, double b,int n)
{

double h=(b-a)/n;
double sum=0.0;
for (int i=0;i<n;i++)
{
sum+=SimpsonSolu(a+i*h, a+(i+1)*h);
}
return sum;
}

/********复化梯形求积分计算函数*************/

double CompositeTrapezoidSolu(double a, double b,int n)
{
double h=(b-a)/n;
double sum=0.0;
for (int i=0;i<n;i++)
{
sum+=TrapezoidSolu(a+i*h, a+(i+1)*h);
}
return sum;
}

/********梯形计算函数*************/
double TrapezoidSolu(double a, double b)
{
return (b-a)/2*(func(a)+func(b));//梯形公式直接求解(误差大)
}

/********辛普生计算函数*************/
double SimpsonSolu(double a, double b)
{
double h=(b-a)/2;
return (b-a)/6*(func(a)+func(b)+4*func(a+h));//辛普生公式直接求解（误差较大）
}

// ConsoleAppMengalotSolu.cpp : 定义控制台应用程序的入口点。
//
/*
*函数功能：函数“x^2”在0到1上的积分
*函数原形：double MengalotSolu(int runs);
*参数：int runs，“打靶”次数
*返回值：积分结果
*时间复杂度：O(n)
*备注：
*日期：2014/11/22
*原创：是
*作者：EbowTang
*Email：tangyibiao520@163.com
*/

#include "stdafx.h"
#include "iostream"
#include <ctime>
#include <iomanip>

using namespace std;

const double pi=3.1415926535897931;
double MengalotSolu(int runs);

int _tmain(int argc, _TCHAR* argv[])
{
cout<<"求函数“x^2”在0到1上的积分"<<endl;
cout<<"/*********随机法求定积分的华丽分割线（打10万次）***********/"<<endl;
int runs=100000;
double res=0.0;
char y='n';
do
{
res= MengalotSolu(runs);
cout<<setprecision(15)<<res;
cout<<"      继续吗？（y/n）";
cin>>y;
} while (y=='y');

return 0;
}

double MengalotSolu(int runs)
{
srand(unsigned int(time(NULL)));//初始化随机函数种子
const int MAX_TIMES = runs;
int inside = 0;
double x=0.0,y=0.0;
for(int i = 0; i < MAX_TIMES; ++i)
{
x =static_cast<double>(rand()) / RAND_MAX;
y =static_cast<double>(rand()) / RAND_MAX;
if(y <= x*x)
++inside;  //统计落入积分区域的点数
}
double res = static_cast<double>(inside) / MAX_TIMES;
return res;
}