1008ac.txt

数值分析的c语言几种算法实现

#include <iostream>
using namespace std;
#define PI 3.1415926535897932 
double Integral(double a, double b, double (*f)(double x, double y, double z),

double eps, double l, double t)
{
    double y[100], q;
    double h ;
	double tot=0.0;
	double hh;  
	h=PI*1.0/t;//步长
	int zz=(int)(1.0*(b-a)/h);
    int num=(int)(h/0.5)+1;
	double inter=h*1.0/(double)num;
	//cout << "zz:" << zz << " num: " << num << " inter: " << inter << endl;
	hh=inter;
	for(int p=1;p<=num;p++)
	{
		hh=inter;
		y[0] = hh *(f(a+(double)(p-1)*inter,l,t) + f(a+(double)p*inter,l,t))/ 2.0;
		int m=1, n=1;
		double ep=eps+1.0;
		while(ep>=eps/(num*zz))
		{
			double pp=0.0;
			for(int i=0;i<n;i++)
			{
				double x = a+(double)(p-1)*inter + (i + 0.5) * hh;
				pp = pp + f(x,l,t);
			}
			pp = (y[0] + hh * pp)/2.0;
			double s=1.0;
			for(int k=1; k<=m; k++)
			{
				s = 4.0 * s;
				q = (s*pp-y[k-1])/(s - 1.0);
				y[k-1] = pp;
				pp = q;
			}
			ep = fabs(q - y[m-1]);
			m = m + 1;
			y[m-1] = q;
			n = n + n;
			hh = hh / 2.0;
		}
		tot+=q;
	}
	//cout << tot << endl;
    tot=tot*zz;
    double sum=tot;
    //cout << sum << " " << tot << endl;
    q=0.0;
	double pp;
    double aa=a+zz*1.0*PI/t;
     hh=b-aa;
     num=(int)(hh/0.3)+1;
	 inter=hh/num;
	//cout << "zz:" << zz << " num: " << num << " inter: " << inter << endl;
	hh=inter;
	for(int p=1;p<=num;p++)
	{
		hh=inter;
		y[0] = hh *(f(aa+(double)(p-1)*inter,l,t) + f(aa+(double)p*inter,l,t))/ 2.0;
		int m=1, n=1;
		double ep=eps+1.0;
		while(ep>=eps/num)
		{
			double pp=0.0;
			for(int i=0;i<n;i++)
			{
				double x = aa+(p-1)*inter + (i + 0.5) * hh;
				pp = pp + f(x,l,t);
			}
			pp = (y[0] + hh * pp)/2.0;
			double s=1.0;
			for(int k=1; k<=m; k++)
			{
				s = 4.0 * s;
				q = (s*pp-y[k-1])/(s - 1.0);
				y[k-1] = pp;
				pp = q;
			}
			ep = fabs(q - y[m-1]);
			m = m + 1;
			y[m-1] = q;
			n = n + n;
			hh = hh / 2.0;
		}
		sum+=q;
	}
	//cout << sum << endl;
    //q*=((b-a)*t*1.0/(PI));
    return(sum/100.0);
}