变步长梯形求解微分方程

变步长梯形求解微分方程2007-03-09 21:16/**
***变步长梯形求解微分方程***

获取初值： 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 算法流程图

**/
#include
#include
#include

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{
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))
printf("final result: %.7f\n",t2); //控制精度
else
{
h=h/2; //修改步长
t1=t2;
goto loop;
}
}
return 0;
}