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

// 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;
}