1. 普通定积分

double f(double x){   //定义被积函数  }

//定义有限n下的的梯形积分，或者辛普森积分

double trap(int n, double a, double b ){      }

//采用变步长技术

double trap_adaptive(double a, double b, double err){

    double s1,s2;

    double n=4;     //initial bin number

    s2 = trap(n,a,b);

    do{            //此处采用do  while语句

        s1=s2;

        n *= 2;

        s2 = trap(n,a,b);

    }while(fabs(s2-s1)>err);

    //return s2;

    return (4*s2-s1)/3;   //combination to improve precision further

}

//主函数调用计算....

2.极限法求解钟摆运动周期

将反常积分归结为：一系列普通定积分的求和。在上述基础代码上，增加对积分范围的选择和对子积分的求和。

以钟摆周期计算题为例子，将积分[0,θ_0]分为:[0,b1], [b1,b2], [b2,b3],.....

其中b_n = θ_0(1-1/2^n)

这部分的核心代码：

double T(double err){

    double s=0,s1,s2;

    double xl=0, xh=theta0/2, dt=0.5;

   s2 = simpson_adaptive(xl, xh, err/10);

    s += s2;

    do{

        s1 = s2;

        xl = xh;

        dt /=2;

        xh = theta0*(1-dt);

        s2 = simpson_adaptive(xl, xh, err/10);

        s += s2;

    }while(fabs(s2-s1)>err);

       return s;

}

习题3中，改变初始角度，计算并输出结果

#include<stdio.h>

#include<math.h>

const double pi=3.14159265358979323846;

double theta0=pi/2;     //在最开始，定义全局变量

......

//在主函数中，改变初始角度，计算并输出结果，用以画图。

int main() {

    double t, err=1e-6;

    FILE *fp = fopen("T.dat","w");

    for(int i=1; i<=8; i++){

        theta0 = pi*pow(2,-i);

        t=T(err);

        fprintf(fp,"%f  %f\n", theta0, t);

    }

    fclose(fp);

    return 0;

}