/**
***变步长梯形求解微分方程***
获取初值： T1 = h/2［ f(x(k)) + f(x(k+1)) ］
n-1
变步长梯形公式： T2n = 1/2*Tn + h/2 * ∑ f ( x(k+1/2) )
k=0
步长： h=b-a/n
属性： 数值积分法
《数值分析简明教程》-2 Editon -高等教育出版社- page 67 算法流程图
代码维护：2005.6.14 DragonLord
**/
#include<iostream.h>
#include<math.h>
#include<stdio.h>
double f(double x)
{
double f;
if(x==0)f=1; // 1
else f=sin(x)/x; // 举例方程 I = ∫ sin(x)/x dx
return f; // 0
}
int main()
{
double a,b,e,h,t1,t2,s,x;
int n=1;
while(cin>>a>>b>>e)
{
h=b-a;
t1=h*(f(a)+f(b))/2; //获取初值
loop: s=0;
x=a+h/2;
while(x<b)
{
s=s+f(x);
x=x+h;
}
t2=t1/2+h*s/2; //求二分后的梯形值
printf("the %d step is %.7f\n",n,t2);
n++;
if((t2-t1<e)&&(t2-t1>-e))
printf("final result: %.7f\n",t2); //控制精度
else
{
h=h/2; //修改步长
t1=t2;
goto loop;
}
}
return 0;
}