/**
变步长Simpson积分
1.获取初值: T1 = h/2[ f(a) + f(b) ],n=1, 步长: h=b-a/n, 且令Sn = Tn
n-1
2.用变步长梯形公式计算: T2n = 1/2*Tn + h/2 * ∑ f ( x(k+1/2) )
k=0
3.用Simpson求积:S2n = (4T2n - Tn ) /3
不满足精度,则加倍分点n,迭代求值.
属性: 数值积分法
《数值计算方法与算法》-2 Editon -科学出版社 P59
《C#数值计算算法编程》-周长发 P315
代码维护:2007.04.20 pengkuny
**/
#include<iostream>
#include<cmath>
![](http://www.cppblog.com/Images/OutliningIndicators/None.gif)
using namespace std;
![](http://www.cppblog.com/Images/OutliningIndicators/None.gif)
#define f(x) (sin(x)) //举例函数
#define epsilon 0.00001 //精度
![](http://www.cppblog.com/Images/OutliningIndicators/None.gif)
//变步长复化梯形公式
double computerAutoT(double aa, double bb)
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedBlockStart.gif)
![](http://www.cppblog.com/Images/OutliningIndicators/ContractedBlock.gif)
{
//迭代初值
long n = 1;
double h = bb-aa; //步长
double t1 = h*(f(aa) + f(bb))/2.0, t2;//t1表示Tn, t2表示T2n
double s1=t1, s2=0; //s1表示Sn, s2表示S2n
double p = epsilon + 1.0;//精度控制
double sum, x;
![](http://www.cppblog.com/Images/OutliningIndicators/InBlock.gif)
while (p >= epsilon)
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{
sum = 0.0;
for (long k=0; k<n; k++)
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{
x = aa + (k+0.5)*h;
sum = sum + f(x);
}
![](http://www.cppblog.com/Images/OutliningIndicators/InBlock.gif)
t2 = (t1 + h*sum)/2.0; //key step
s2 = (4.0*t2 - t1)/3.0; //key step
p = fabs(s2-s1);
t1 = t2; s1 = s2; n = n+n; h = h/2.0;
}
![](http://www.cppblog.com/Images/OutliningIndicators/InBlock.gif)
cout<<"最终分点n:"<<n<<endl;
return (s2);
}
![](http://www.cppblog.com/Images/OutliningIndicators/None.gif)
![](http://www.cppblog.com/Images/OutliningIndicators/None.gif)
int main()
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedBlockStart.gif)
![](http://www.cppblog.com/Images/OutliningIndicators/ContractedBlock.gif)
{
double a,b;
cout<<"变步长复化梯形积分,请输入积分范围a,b:"<<endl;
cin>>a>>b;
![](http://www.cppblog.com/Images/OutliningIndicators/InBlock.gif)
cout<<"积分结果:"<<computerAutoT(a, b)<<endl;
![](http://www.cppblog.com/Images/OutliningIndicators/InBlock.gif)
system("pause");
return 0;
}
posted on 2007-04-20 10:55
哈哈 阅读(1130)
评论(0) 编辑 收藏 引用