关于数值方法的一些算法解析（3）

再次继续~最后了

1、复合求积公式计算定积分

//

 

#include <iostream>

#include <iomanip>

#include <cmath>

using namespace std;

//

#define e 2.718281828459

#define epsilon 0.0001  //精度

#define MAXREPT  10   //迭代次数,到最后仍达不到精度要求,则输出T(m=10).

void Print() 

{

system("cls");

cout<<"/t0) 退出!"<<endl;

cout<<"/t1) 复合梯形求解!"<<endl;

cout<<"/t2) 复合抛物线求解!"<<endl;

cout<<"/t3) 龙贝格求解!"<<endl;

cout<<"/t你的选择：";

}

double f(double x)

{

return x*pow(e,x);

}

//举例函数

using namespace std;

int main()

{

void Compos_Trapezoidal();

void Compo_Parabolic();

void Romberg();

int Choose;

bool flag=true;

while(flag)

{

Print();

cin>>Choose;

if(Choose!=0&&Choose!=1&&Choose!=2&&Choose!=3){

cout<<"/tError: unvalid choose!Please Choose again!"<<endl;

system("pause");

continue;

}

//

switch(Choose)

{

case 0:flag=false;break;

case 1:Compos_Trapezoidal();break;

case 2:Compo_Parabolic();break;

case 3:Romberg();break;

}

}

return 0;

}

///

void Compos_Trapezoidal()

{

int k,n;

double a,b,h;

double Sum=0.0;

cout<<"请输入区间端点a和b：";

cin>>a>>b;

cout<<"请输入等分数(n):";

cin>>n;

h=(b-a)/n;

for(k=1;k<n;k++){

//

Sum+=f(a+h*k);

}

Sum*=h;

Sum+=h*(f(a)+f(b))/2;

cout<<"复合梯形求解为:"<<setprecision(8)<<fixed<<Sum<<endl;

cout<<"e的次方结果为:"<<pow(e,2)<<endl;

cout<<"两者的误差为:"<<pow(e,2)-Sum<<endl;

system("pause");

}

//

void Compo_Parabolic()

{

int n;

int i;

double a,b,h;

double Sum=0.0;

double Temp_1(0),Temp_2(0);

cout<<"请输入区间端点a和b：";

cin>>a>>b;

cout<<"请输入等分数(n):";

cin>>n;

h=(b-a)/(2*n);

//

for(i=0;i<n;i++){

Temp_1+=f(a+(2*i+2)*h);

Temp_2+=f(a+(2*i)*h);

}

Sum+=4*Temp_1+2*Temp_2;

cout<<Sum<<endl;

Sum+=f(a)+f(b);

Sum*=h/3;

cout<<"复合抛物线求解为:"<<setprecision(8)<<fixed<<Sum<<endl;

cout<<"e的次方结果为:"<<pow(e,2)<<endl;

cout<<"两者的误差为:"<<pow(e,2)-Sum<<endl;

system("pause");

}

//

void Romberg()

{

int m, n;//m控制迭代次数, 而n控制复化梯形积分的分点数. n=2^m

    double h, x;

    double s, q;

    double ep; //精度要求

    double *y = new double[MAXREPT];//为节省空间,只需一维数组

    //每次循环依次存储Romberg计算表的每行元素,以供计算下一行,算完后更新

    double p;//p总是指示待计算元素的前一个元素(同一行)

double aa,bb;

cout<<"请输入区间端点a和b：";

cin>>aa>>bb;

    //迭代初值

    h = bb - aa;

    y[0] = h*(f(aa) + f(bb))/2.0;

    m = 1;

    n = 1;

    ep = epsilon + 1.0;

 

    //迭代计算

    while ((ep >= epsilon) && (m < MAXREPT))

    {

        //复化积分公式求T2n(Romberg计算表中的第一列),n初始为1,以后倍增

        p = 0.0;

        for (int i=0; i<n; i++)//求Hn

        {

            x = aa + (i+0.5)*h;

            p = p + f(x);

        }

        p = (y[0] + h*p)/2.0;//求T2n = 1/2(Tn+Hn),用p指示

 

        //求第m行元素,根据Romberg计算表本行的前一个元素(p指示),

        //和上一行左上角元素(y[k-1]指示)求得.        

        s = 1.0;

        for (int k=1; k<=m; k++)

        {

            s = 4.0*s;

            q = (s*p - y[k-1])/(s - 1.0);

            y[k-1] = p;

            p = q;

        }

 

        p = fabs(q - y[m-1]);

        m = m + 1;

        y[m-1] = q;

        n = n + n; h = h/2.0;

    }

cout<<"龙贝格求解为:"<<setprecision(8)<<fixed<<q<<endl;

cout<<"e的次方结果为:"<<pow(e,2)<<endl;

cout<<"两者的误差为:"<<pow(e,2)-q<<endl;

system("pause");

}