关于数值方法的一些算法解析（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");
}