数值计算数值积分实现

原创作品,出自 “晓风残月xj” 博客,欢迎转载,转载时请务必注明出处(http://blog.csdn.net/xiaofengcanyuexj)。

由于各种原因,可能存在诸多不足,欢迎斧正!

一、分工合作

    本次报告组长徐进负责写程序,小组成员张映红、邓洋、潘传军、杨林岱、王一鹤提出意见与修改建议,邓洋负责写报告word文档,杨林岱负责做PPT,黄江、易健、连王强等着手word文档与PPT的后期制作,最后由资深发言人王一鹤负责组间展示与交流。

二、引言

     在大量的实际问题中,定积分的计算是需要的。由微积分基本定理可知,如果能找到被积函数fx)的原函数Fx), 则用牛顿-莱布尼兹公式可得

                 I[f]=fab(x)dx=F(b)-F(a).

但是牛顿-莱布尼兹公式并不能解决所有函数积分问题,因为大量的被积分函数,如sinx/xe-x2等,无法用初等函数表示其原函数的,此外被积函数为一组数据时,牛顿-莱布尼兹公式更是无能为力。但从理论上看,定积分为一确定的数值是客观存在的,可以尝试用计算机依据定义进行数值积分。当然我们讨论的都是定积分,定积分的基本思想为大化小,常代变,近似和,取极限;我们知道,计算机处理的是离散的、不连续的问题,这与数值积分的思想高度吻合,因此有必要运用计算机研究定积分的数值计算问题。

三、算法描述(Linear equation

1Constant_rule

使用数组nun_interval[j]来保存将积分区间分的段数2^j,用整形变量inc保存每一段小区间长度,用num_f[j]来记录记录对应函数值的个数,并且此时num_f[j]=nun_interval[j]+1。然后根据choice的值来选择是进行课本函数积分的验证,还是通过输入函数进行识别。如果是识别用户输入函数,则需要使用函数post_Exp()将用户输入函数转换为编译器可以识别的代码,并执行sum[j]+=inc*post_Exp(x0+i*inc);否则,根据Constant公式:

注释:(此公式是在将积分区间n等分之后,取每一个小块,将其看成是矩形,利用矩形面积公式S=长(该小段左端点处的函数值)*宽(小段的长度)求每个小矩形的面积并进行累加来完成积分运算的,并且取得的对应函数的点的个数总是比积分区间的段数多1,当n的值越大,积分结果距离真实值的误差越小,当n值足够大时,积分值越来越趋近真实值)

直接执行sum[j]+=inc*sin(x0+i*inc)的代码段,其中sum[j]用来保存当前的积分值。然后调用函数print()打印出:积分区间分段个数、积分误差、对应函数值的个数、积分结果等数据。

2、Trapezium_rule

    使用数组nun_interval[j]来保存将积分区间分的段数2^j,用整形变量inc保存每一段小区间长度,用num_f[j]来记录记录对应函数值的个数,并且此时num_f[j]=nun_interval[j]+1。然后根据choice的值来选择是进行课本函数积分的验证,还是通过输入函数进行识别。如果是识别用户输入函数,则需要使用函数post_Exp()将用户输入函数转换为编译器可以识别的代码,并执行sum[j]+=inc*(post_Exp(x0+i*inc)+post_Exp(x0+(i+1)*inc))/2.0;否则,根据Trapezium公式:

注释:(此公式是在将积分区间n等分之后,取每一个小块,将其看成是梯形,利用梯形面积公式S=(上底(该小段左端点处的函数值)+下底(该小段右端点处的函数值))*高(小段的长度)/2求每个小梯形的面积并进行累加来完成积分运算的,并且取得的对应函数的点的个数总是比积分区间的段数多1,当n的值越大,积分结果距离真实值的误差越小,当n值足够大时,积分值越来越趋近真实值)

直接执行sum[j]+=inc*(sin(x0+i*inc)+sin(x0+(i+1)*inc))/2.0;的代码段,其中sum[j]用来保存当前的积分值。然后调用函数print()打印出:积分区间分段个数、积分误差、对应函数值的个数、积分结果等数据。

3、Mid_point_rule

    使用数组nun_interval[j]来保存将积分区间分的段数2^j,用整形变量inc保存每一段小区间长度,用num_f[j]来记录记录对应函数值的个数,并且此时num_f[j]=nun_interval[j]。然后根据choice的值来选择是进行课本函数积分的验证,还是通过输入函数进行识别。如果是识别用户输入函数,则需要使用函数post_Exp()将用户输入函数转换为编译器可以识别的代码,并执行sum[j]+=inc* (post_Exp( (x0+i*inc+x0+(i+1)*inc) /2.0 ));否则,根据Mid_point公式:

     注释:(此公式是在将积分区间n等分之后,取每一个小块,将其看成是矩形,利用矩形面积公式S=长(该小段中点处的函数值)*宽(小段的长度)求每个小矩形的面积并进行累加来完成积分运算的,并且取得的对应函数的点的个数总是比积分区间的段数是一样的·,当n的值越大,积分结果距离真实值的误差越小,当n值足够大时,积分值越来越趋近真实值)

直接执行sum[j]+=inc* (sin( (x0+i*inc+x0+(i+1)*inc) /2.0 ))的代码段,其中sum[j]用来保存当前的积分值。然后调用函数print()打印出:积分区间分段个数、积分误差、对应函数值的个数、积分结果等数据。

4、Simpson_rule

    使用数组nun_interval[j]来保存将积分区间分的段数2^j,用整形变量inc保存每一段小区间长度,用num_f[j]来记录记录对应函数值的个数,并且此时num_f[j]=nun_interval[j]*2+1。然后根据choice的值来选择是进行课本函数积分的验证,还是通过输入函数进行识别。如果是识别用户输入函数,则需要使用函数post_Exp()将用户输入函数转换为编译器可以识别的代码,并执行sum[j-1]+=inc/3.0 *(post_Exp(x0+i*inc)+4*post_Exp(x0+(i+1)*inc)+post_Exp(x0+(i+2)*inc));否则,根据Simpson公式:

    注释:(此公式是在将积分区间n等分之后,取每一个小段,每次取每一个等分点的函数值以及接下来一个等分点的函数值的4倍,还有在下一个等分点的函数值的平均值,与每小段的长度相乘得到每一小块的面积,并累加求得整个积分区间上的积分值,并且取得的对应函数的点的个数总是比积分区间的段数的2倍加1,当n的值越大,积分结果距离真实值的误差越小,当n值足够大时,积分值越来越趋近真实值)

直接执行sum[j-1]+=inc/3.0 *(sin(x0+i*inc)+4*sin(x0+(i+1)*inc)+sin(x0+(i+2)*inc));的代码段,其中sum[j]用来保存当前的积分值。然后调用函数print()打印出:积分区间分段个数、积分误差、对应函数值的个数、积分结果等数据。

5Gauss_quadrature

    使用数组nun_interval[j]来保存将积分区间分的段数2^j,用整形变量inc保存每一段小区间长度,用 num_f[j]来记录记录对应函数值的个数,并且此时num_f[j]=nun_interval[j]*2。然后根据choice的值来选择是进行课本函数积分的验证,还是通过输入函数进行识别。如果是识别用户输入函数,则需要使用函数post_Exp()将用户输入函数转换为编译器可以识别的代码,并执行sum[j]+=inc/2.0 * ( post_Exp(x0+inc*i+(1.0/2.0-sqrt(3.0)/6.0)*inc) + post_Exp(x0+inc*i+(1.0/2.0+sqrt(3.0)/6.0)*inc) );否则,根据Gauss_quadrature公式:

    注释:(此公式是在将积分区间n等分之后,取每一个小段,取比每一个等分点小(1/2-/6以及比这个等分点大(1/2+/6处的函数值的平均值,乘以每小段的长度为每一小块的面积,并累加求得整个积分区间上的积分值,并且取得的对应函数的点的个数总是比积分区间的段数的2倍,当n的值越大,积分结果距离真实值的误差越小,当n值足够大时,积分值越来越趋近真实值)

直接执行sum[j]+=inc/2.0 * ( sin(x0+inc*i+(1.0/2.0-sqrt(3.0)/6.0)*inc) + sin(x0+inc*i+(1.0/2.0+sqrt(3.0)/6.0)*inc) );的代码段,其中sum[j]用来保存当前的积分值。然后调用函数print()打印出:积分区间分段个数、积分误差、对应函数值的个数、积分结果等数据。

6、Gauss_quadrature_three_point

    使用数组nun_interval[j]来保存将积分区间分的段数2^j,用整形变量inc保存每一段小区间长度,用num_f[j]来记录记录对应函数值的个数,并且此时num_f[j]=nun_interval[j]*3。然后根据choice的值来选择是进行课本函数积分的验证,还是通过输入函数进行识别。如果是识别用户输入函数,则需要使用函数post_Exp()将用户输入函数转换为编译器可以识别的代码,并执行sum[j-1]+=inc/9.0 * ( 8*post_Exp(x0+inc+inc*i) + 5*post_Exp(x0+inc+inc*i-sqrt(3.0/5.0)*inc)  +5*post_Exp(x0+inc+inc*i+sqrt(3.0/5.0)*inc) );否则,根据Gauss_quadrature_three_point公式:

     注释:(此公式是在将积分区间n等分之后,取每一个小段,每个小块的面积需要求到三个点的函数值,分别是x0点的,点的,以及点的,每一小块的面积如公式给出的方法求出,并累加求得整个积分区间上的积分值,并且取得的对应函数的点的个数总是比积分区间的段数的3倍,当n的值越大,积分结果距离真实值的误差越小,当n值足够大时,积分值越来越趋近真实值)

直接执行sum[j-1]+=inc/9.0 * ( 8*sin(x0+inc+inc*i) + 5*sin(x0+inc+inc*i-sqrt(3.0/5.0)*inc)  +5*sin(x0+inc+inc*i+sqrt(3.0/5.0)*inc))的代码段,其中sum[j]用来保存当前的积分值。然后调用函数print()打印出:积分区间分段个数、积分误差、对应函数值的个数、积分结果等数据。

四、算法分析与运行结果分析

  1.对于Constant_rule法,每段小区间用矩形表示,比较粗糙,误差为O^2),系数为f的一阶倒数。

  

从运算结果收敛速度较慢

      

    2.对于Trapezium_rule法,每段小区间用梯形表示,一般而言较Constant_rule有所改善,比较粗糙,误差为O^2),系数为f的二阶倒数。

     3.对于Mid_point_rule法,每段小区间用中点函数值表示,有所改善,比较粗糙,误差为O^2),系数为f的二阶倒数。

 

     4.对于Simpson_rule法,每段小区间用三个函数值表示,误差为O(^4),收敛速度也比较快。

    

5.对于Gauss_quadrature法,每段小区间用三个函数值表示,

误差为O(^3,收敛速度也比较快。

   

6.对于Gauss_quadrature_three_point法,每段小区间用三个函数值表示,误差为O(^6),而收敛速度最快。

     

五、特色

     1.实现了用户输入函数

      我们的程序可以根据用户的具体计算需要输入函数(由于写程序初期有些问题未意识到,所以并不能覆盖所有的函数,但较之在代码中封装单一sinx)函数有很大改善),然后通过自己写的头文件"识别函数程序3.h",完成函数的识别工作。具体工作为以字符串形式输入用户函数,然后结合表达式的转换与计算知识将输入的中缀表达式转换为后缀表达式存储在自定义堆栈中,在此基础上根据用户输入的积分区间,调用<math>库中的函数计算后缀表达式的值(具体实现见源代码)。

   我们组可以识别的函数举例(有些函数可能不常见,这里不做解释)

fabs     cos     sin      tan      asin       acos     atan  

sinh    cosh     tanh    exp      floor       log     

以及它们与x的任意组和,复杂多项式等等,在注意了输入格式的前提下基本可以实现绝大多数函数的数值积分。

     下面是一组用户随机函数运用课本数值积分的6中思想得到的结果,有兴趣的同学可以验证一下:

    用户手动输入随机函数:absin(x)^cos(x)-x*(1-x^2+10*x*x^2*cos(x))-loge(x)

有各种函数与多项式,其中积分区间用户输入,此处a=2,b=3

  

   6种积分方法的收敛速度不同,但数值积分结果逼近369.371,相互验证了各自结果的正确性,从而说明我们组数值积分方法的实现程序和函数识别程序没有明显错误。

     2.菜单界面

        本程序设计了菜单界面,可以选择使用123号功能,1功能代表直接对absinx=2积分结果的验证,2功能则是对用户输入的函数进行求积。之所以要分开,主要原因有两点:

    A\1号功能absinx=2是在知道积分结果的情况下的验证性计算,需输出误差与收敛速度,而2号功能是自己随机输入函数的求解性计算,在不知道真实值的前提下是无法直到误差与收敛速度的;

B\1号功能直接调用<math>库中的sin(x)函数,资源占用小;而2号功能-自己输入的函数则要通过繁琐的字符串识别与表达式的计算,涉及进出栈需占用大量系统资源,时间和空间尤其是时间资源消耗较大,有一定运行时间,所以分开便于展示。

关于菜单,我们积极响应曾老师的号召运用英语做的,

虽然有点班门弄斧的意味,但不失为一次大胆学以致用运用英语的尝试。

六、总结

   数值积分在科学和工程中占有很高的地位,方法也很多,是解决定积分的重要方法。Constant ruleTrapezium ruleMid point ruleSimpson ruleGauss quadratureGauss quadrature three point等都是利用定积分的基本思想:大化小,常代变,近似和,取极限;结合计算机处理的问题是离散的、不连续的特点,用计算机有效的处理数值积分问题。Constant ruleTrapezium ruleMid point ruleSimpson ruleGauss quadratureGauss quadrature three point等规则为数值积分提供了理论依据(包括误差,收敛速度等),而计算机则为数值积分提供了实现工具。这是理论与技术的完美结合。

     通过这次以及前几次的计算方法的学习与交流,我们加深了对计算机解决实际问题方法的认识:

     1.提出问题

     2.建立模型

     3.理论分析(这点很重要,就拿这次而言,这6种数值积分的基本思想,为数值积分的实现提供了理论依据,同时也证明了其方法的正确性与可行性,使我们编程时有了明确的指导思想不至于盲干)

     4.编程实现


//头文件
#include<iostream>
#include<math.h>
#include<string>
#include<stack>
using namespace std;

char expression[200];

bool operator_priority(char temp_operator,char expression)
//若栈顶运算符temp_operator优先级比当前表达式的高,则返回true,否则返回false
{
    if(temp_operator=='-'||temp_operator=='+')
    {
        if(expression=='-'||expression=='+'||expression==')'||expression=='#')
            return true;
        else if(expression=='/'||expression=='*'||expression=='('||expression=='^')
            return false;
    }
    
    else if(temp_operator=='/'||temp_operator=='*')
    {
        if(expression=='-'||expression=='+'||expression=='/'||expression=='*'||expression==')'||expression=='#')
            return true;
        else if(expression=='('||expression=='^')
            return false;
    }
    
    else if(temp_operator=='(')
    {
        if(expression!=')')
            return false;
		else true;
    }
    
    else if(temp_operator==')')
        return true;
    
    else if(temp_operator=='#')
        return false;
	
	else if(temp_operator=='^')/
	{
		if(expression=='(')
			return false;
     	else true;
	}
}

bool is_operator(char expression)
//判断expression是否为运算符,若是返回true,否则返回false
{
	if(expression=='-'||expression=='+'||expression=='*'||expression=='/'||expression==')'||expression=='('||expression=='^'||expression=='#')
		return true;
	else return false;
}

int len;
string post;
string str[200];
void mid_trans_post()
//将中缀表达式转化为后缀表达式存储在str[200]
{
    int i;
    gets(expression);
    stack<char>Stack_operator;
    Stack_operator.push('#');
    i=0;
	//  cout<<"expression="<<expression<<endl;
    
    
    while(i<strlen(expression))
    {
        if(expression[i]>='0'&&expression[i]<='9')//代表的是数据
            post=post+expression[i];
        
        else if(expression[i]>='a'&&expression[i]<='z')//代表的是数据//
            post=post+expression[i];
        
		else if(i>0&&!is_operator(expression[i-1])&&is_operator(expression[i])){str[len++]=post;/*cout<<" post="<<post;*/post.erase(0,post.size());}
				
		
		if(is_operator(expression[i]))//代表是操作符
        {
leap:   char temp_operator=Stack_operator.top();        
        if(temp_operator=='('&&expression[i]==')')
        {   
            Stack_operator.pop();
        }
        else if(temp_operator=='#'&&expression[i]=='#')
        {Stack_operator.pop();break;}
        else if(temp_operator==')'&&expression[i]=='(')
        {   cout<<"你的表达是有错!"<<endl;exit(0);}
        else if(!(operator_priority(temp_operator,expression[i])))
        {
            Stack_operator.push(expression[i]);
        }
        else if(operator_priority(temp_operator,expression[i]))
        {   
            Stack_operator.pop();
            str[len++]=temp_operator;
            while(!Stack_operator.empty())
            {
                goto leap;
            }
            
        }
        }
        i++;
    }    
}


double post_Exp(double iteration_x)
//计算代值iteration_x之后后缀表达式的值
{
    double x1,x,x2;
    int i;
    stack<double>Stack_operand;
    
    
    for(i=0;i<len;i++)
    {
        x=0;
        if(str[i][0]>='0'&&str[i][0]<='9')
        {
            int j=0;
            while(j<str[i].size())
            {
                x=x*10+str[i][j]-'0';
                j++;
            }
			// cout<<" x="<<x<<" ";
            Stack_operand.push(x);
        }
        else if(str[i][0]>='a'&&str[i][0]<='z')
        {
			double temp;
			if(str[i]=="x")temp=iteration_x;
			else if(str[i]=="fabs")temp=fabs(iteration_x);
			else if(str[i]=="cos")temp=cos(iteration_x);
            else if(str[i]=="sin")temp=sin(iteration_x);
            else if(str[i]=="tan")temp=tan(iteration_x);
            else if(str[i]=="asin")temp=asin(iteration_x);
            else if(str[i]=="acos")temp=acos(iteration_x);
            else if(str[i]=="atan")temp=atan(iteration_x);
			//  else if(str[i]=="atan2")temp=atan2(1,2);
            else if(str[i]=="sinh")temp=sinh(iteration_x);
            else if(str[i]=="cosh")temp=cosh(iteration_x);
            else if(str[i]=="tanh")temp=tanh(iteration_x);
            else if(str[i]=="exp")temp=exp(iteration_x);
            else if(str[i]=="floor")temp=floor(iteration_x);
			//  else if(str[i]=="fmod")temp=fmod(5,3);
            else if(str[i]=="log")temp=log(iteration_x);
			//  else if(str[i]=="log10")temp=log10(105);
			Stack_operand.push(temp);
        }
        else 
        {
            x2=Stack_operand.top();
			//   cout<<" x2="<<x2<<" ";
            Stack_operand.pop();
			x1=Stack_operand.top();
			//  cout<<" x1="<<x1<<" ";
            Stack_operand.pop();
            
            switch(str[i][0])
            {
                
            case '+':x1+=x2;break;
            case '-':x1-=x2;break;
            case '*':x1*=x2;break;
            case '/':x1/=x2;break;
			case '^':x1=pow(x1,x2);
            }
            Stack_operand.push(x1);
        }
    }
    return Stack_operand.top();    
}


void menu()
//输出提示并调用相关函数进行计算
{
	int i;
	cout<<"*****************************************************************"<<endl;
	cout<<"Please input the  function(care for the format):"<<endl;
	
	len=0;
	mid_trans_post();
	cout<<"followings is the infix of the  function:"<<endl;
	for(i=0;i<strlen(expression)-1;i++)
	{
		if(expression[i]!=' ')cout<<expression[i];
	}
		
		cout<<endl<<"followings is the suffix of the  function:"<<endl;
		for(i=0;i<len;i++)
		{
			cout<<str[i]<<" ";
		}
	cout<<endl<<"*****************************************************************"<<endl;
		
}
/*
int main()
{
    while(1)
    {
		menu();
		
		getchar();
    }
    
    return 0;
}*/

//源文件
#include<iostream>
#include"识别函数程序3.h"
#include<iomanip>
#include<math.h>
#define MAX 100
using namespace std;

//测试实例
//      sin(x)*dx|0——PI=2

#define PI 3.141592653
double x0=0;
double x1=PI;/赋初值

double inc;/表示小区间长度
double sum[MAX];表示积分结果
double e[MAX];//记录积分误差
int nun_interval[MAX];记录区间个数
int num_f[MAX];//记录对应函数值的个数
int choice;//表示是进行课本函数积分的验证,还是通过输入函数进行识别

void print(int num)
//以固定格式输出
{
	int i;

	/输出测试实例的结果,带误差的比较
	if(1==choice)
	{
		for(i=0;i<num;i++)
		{
			e[i]=sum[i]-2;
		}
		for(i=0;i<num-1;i++)
		{
			cout<<setw(15)<<nun_interval[i]<<setw(15)<<num_f[i]<<setw(20)<<sum[i]-2.0<<setw(20)<<e[i+1]/e[i]<<endl;
		}
		cout<<setw(15)<<nun_interval[i]<<setw(15)<<num_f[i]<<setw(20)<<sum[i]-2.0<<setw(20)<<endl;
	}

	/输出输出随机函数的结果,由于不知道精确值,所以不带误差的比较
    else if(2==choice)
	{
		for(i=0;i<num-1;i++)
		{
			cout<<setw(15)<<nun_interval[i]<<setw(15)<<num_f[i]<<setw(20)<<sum[i]<<endl;
		}
		cout<<setw(15)<<nun_interval[i]<<setw(15)<<num_f[i]<<setw(20)<<sum[i]<<endl;
	}
}

double Constant_rule()
//Constant rule确定小区间面积
{
	int j;
	if(1==choice)
		cout<<setw(15)<<"No.interval"<<setw(15)<<"No.f(x)"<<setw(20)<<"Constant rule"<<setw(20)<<"e[2*n]/e[n]"<<endl;
	else cout<<setw(15)<<"No.interval"<<setw(15)<<"No.f(x)"<<setw(20)<<"Result of Constant rule "<<endl;

	///反复进行,逐步将区间长度减小,逐步逼近精确值
	for(j=0;j<20;j++)
	{
		nun_interval[j]=(int)pow(2,j);
		sum[j]=0;
		inc=(x1-x0)/nun_interval[j];
		int i=0;

		利用递推法求当前区间分法的积分值sum[j]
		while(i<nun_interval[j])
		{
			if(2==choice)sum[j]+=inc*post_Exp(x0+i*inc);
			else sum[j]+=inc*sin(x0+i*inc);
			i++;
		}
		num_f[j]=nun_interval[j]+1;
	}
	print(j);
	return sum[j-1];
}

double Trapezium_rule()
//Trapezium rule确定小区间面积
{
	int j;
	if(1==choice)
		cout<<setw(15)<<"No.interval"<<setw(15)<<"No.f(x)"<<setw(20)<<"Trapezium rule"<<setw(20)<<"e[2*n]/e[n]"<<endl;
	else cout<<setw(15)<<"No.interval"<<setw(15)<<"No.f(x)"<<setw(20)<<"Result of Trapezium rule "<<endl;

	///反复进行,逐步将区间长度减小,逐步逼近精确值
	for(j=0;j<20;j++)
	{
		nun_interval[j]=(int)pow(2,j);
		sum[j]=0;
		inc=(x1-x0)/nun_interval[j];
		int i=0;

		利用递推法求当前区间分法的积分值sum[j]
		while(i<nun_interval[j])
		{
			if(2==choice)sum[j]+=inc*(post_Exp(x0+i*inc)+post_Exp(x0+(i+1)*inc))/2.0;
			else sum[j]+=inc*(sin(x0+i*inc)+sin(x0+(i+1)*inc))/2.0;
			i++;
		}
		num_f[j]=nun_interval[j]+1;
	}
	print(j);
	return sum[j-1];	
}

double Mid_point_rule()
//Mid_point rule确定小区间面积
{
	int j;
	if(1==choice)
		cout<<setw(15)<<"No.interval"<<setw(15)<<"No.f(x)"<<setw(20)<<"Mid_point rule"<<setw(20)<<"e[2*n]/e[n]"<<endl;
	else cout<<setw(15)<<"No.interval"<<setw(15)<<"No.f(x)"<<setw(20)<<"Result of Mid_point rule "<<endl;

	///反复进行,逐步将区间长度减小,逐步逼近精确值
	for(j=0;j<20;j++)
	{
		nun_interval[j]=(int)pow(2,j);
		sum[j]=0;
		inc=(x1-x0)/nun_interval[j];
		int i=0;

		利用递推法求当前区间分法的积分值sum[j]
		while(i<nun_interval[j])
		{
			if(2==choice)sum[j]+=inc* (post_Exp( (x0+i*inc+x0+(i+1)*inc) /2.0 ));
			else sum[j]+=inc* (sin( (x0+i*inc+x0+(i+1)*inc) /2.0 ));
			i++;
		}
		num_f[j]=nun_interval[j];
	}
	print(j);
	return sum[j-1];
}

double Simpson_rule()
//Simpson's_rule()确定小区间面积
{
	int j;
	if(1==choice)
		cout<<setw(15)<<"No.interval"<<setw(15)<<"No.f(x)"<<setw(20)<<"Simpson's rule"<<setw(20)<<"e[2*n]/e[n]"<<endl;
	else cout<<setw(15)<<"No.interval"<<setw(15)<<"No.f(x)"<<setw(20)<<"Result of Simpson's rule "<<endl;

	///反复进行,逐步将区间长度减小,逐步逼近精确值
	for(j=0;j<11;j++)
	{
		nun_interval[j]=(int)pow(2,j);
		sum[j]=0;
		inc=(x1-x0)/nun_interval[j];
		int i=0;

		利用递推法求当前区间分法的积分值sum[j]
		while(j!=0&i<nun_interval[j])
		{
			if(2==choice)sum[j-1]+=inc/3.0 *(post_Exp(x0+i*inc)+4*post_Exp(x0+(i+1)*inc)+post_Exp(x0+(i+2)*inc));
			else sum[j-1]+=inc/3.0 *(sin(x0+i*inc)+4*sin(x0+(i+1)*inc)+sin(x0+(i+2)*inc));
			i=i+2;
		}
		num_f[j]=nun_interval[j]*2+1;
	}
	print(j-1);
	return sum[j-2];
}


double Gauss_quadrature()
//Gauss_quadrature()确定小区间面积
{
	int j;
	if(1==choice)
		cout<<setw(15)<<"No.interval"<<setw(15)<<"No.f(x)"<<setw(20)<<"Gauss quadrature"<<setw(20)<<"e[2*n]/e[n]"<<endl;
	else cout<<setw(15)<<"No.interval"<<setw(15)<<"No.f(x)"<<setw(20)<<"Result of Gauss quadrature "<<endl;

	///反复进行,逐步将区间长度减小,逐步逼近精确值
	for(j=0;j<10;j++)
	{
		nun_interval[j]=(int)pow(2,j);
		sum[j]=0;
		inc=(x1-x0)/nun_interval[j];
		int i=0;

		利用递推法求当前区间分法的积分值sum[j]
		while(i<nun_interval[j])
		{
			if(2==choice)sum[j]+=inc/2.0 * ( post_Exp(x0+inc*i+(1.0/2.0-sqrt(3.0)/6.0)*inc) + post_Exp(x0+inc*i+(1.0/2.0+sqrt(3.0)/6.0)*inc) );
			else sum[j]+=inc/2.0 * ( sin(x0+inc*i+(1.0/2.0-sqrt(3.0)/6.0)*inc) + sin(x0+inc*i+(1.0/2.0+sqrt(3.0)/6.0)*inc) );
			i++;
		}
		num_f[j]=nun_interval[j]*2;
	}
	print(j);
	return sum[j-1];
}




double Gauss_quadrature_three_point()
//Gauss_quadrature_three_point()确定小区间面积
{
	int j;
	if(1==choice)
		cout<<setw(15)<<"No.interval"<<setw(15)<<"No.f(x)"<<setw(20)<<"Gauss quadrature three point"<<setw(20)<<"e[2*n]/e[n]"<<endl;
	else cout<<setw(15)<<"No.interval"<<setw(15)<<"No.f(x)"<<setw(20)<<"Result of Gauss quadrature three point "<<endl;

	///反复进行,逐步将区间长度减小,逐步逼近精确值
	for(j=0;j<10;j++)
	{
		nun_interval[j]=(int)pow(2,j);
		sum[j]=0;
		inc=(x1-x0)/nun_interval[j];                                                                                                                                                                                                                                                                                                                          
		num_f[j]=3*nun_interval[j];
		int i=0;  
		
		利用递推法求当前区间分法的积分值sum[j]
		while(j!=0&&i<nun_interval[j])
		{
			if(2==choice)sum[j-1]+=inc/9.0 * ( 8*post_Exp(x0+inc+inc*i) + 5*post_Exp(x0+inc+inc*i-sqrt(3.0/5.0)*inc)  +5*post_Exp(x0+inc+inc*i+sqrt(3.0/5.0)*inc) );
			else sum[j-1]+=inc/9.0 * ( 8*sin(x0+inc+inc*i) + 5*sin(x0+inc+inc*i-sqrt(3.0/5.0)*inc)  +5*sin(x0+inc+inc*i+sqrt(3.0/5.0)*inc) );
			i=i+2;
		}
		//sum[j]/=2;
	}8
	print(j-1);
	return sum[j-2];
}



int main()
{
	
	while(1)
	{
		/根据菜单提示进行相应操作///
		int flag=0;
		do
		{
			if(flag)cout<<"Your choice is illegal,please choose again!"<<endl;
			cout<<"============================menu================================"<<endl;
			cout<<"              1.Integral verification on the books"<<endl;
			cout<<"              2.Input new function and do numerical integration"<<endl;
			cout<<"              3.Quits the program"<<endl;
			cout<<"Please choose [ ]\b\b";
			cin>>choice;
			getchar();
			cout<<"================================================================"<<endl<<endl;
			flag++;
		}while(choice!=1&&choice!=2&&choice!=3);	
		if(2==choice)	
		{
			menu();
			cout<<endl<<"Please input interval [x0,x1]:";
			cin>>x0>>x1;
			cout<<"followings is the root of the  function:"<<endl<<endl;
		}
		else if(3==choice)
			exit(0);
		
		double result;
        cout<<"*********here is the result of Constant rule***********"<<endl;
		result=Constant_rule();
		cout<<"result of Constant rule:"<<result<<endl<<endl;

		cout<<"*********here is the result of Trapezium rule***********"<<endl;
		result=Trapezium_rule();
		cout<<"result of Trapezium rule:"<<result<<endl<<endl;

		cout<<"*********here is the result of Mid point rule***********"<<endl;
		result=Mid_point_rule();
		cout<<"result of Mid point rule:"<<result<<endl<<endl;

		cout<<"*********here is the result of Simpson rule***********"<<endl;
		result=Simpson_rule();
		cout<<"result of Simpson rule:"<<result<<endl<<endl;

		cout<<"*********here is the result of Gauss quadrature*********"<<endl;
		result=Gauss_quadrature();
		cout<<"result of Gauss quadrature:"<<result<<endl<<endl;

		cout<<"*********here is the result of Gauss quadrature three point**********"<<endl;
		result=Gauss_quadrature_three_point();
		cout<<"result of Gauss quadrature three point:"<<result<<endl<<endl;
	}
	return 0;	
}


   以上程序有很多不足之处,望多包涵!欢迎交流!

评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值