梯形多步法和辛普森积分
具体算法以及相关定理参考以下链接:
#include <stdio.h>
#include <string.h>#include <stdlib.h>
#include <math.h>
#define EPS 1.0E-14 //计算精度
#define DIV 1000 //分割区间,值越大,精确值越高
#define ITERATE 20 //二分区间迭代次数,区间分割为2^n,初始化应该小一点,否则会溢出
#define RECTANGLE 0 //矩形近似
#define TRAPEZOID 1 //梯形近似
#define TRAPEZOID_FORMULA 2 //递推梯形公式
#define SIMPSON_FORMULA 3 //辛普森公式
#define BOOL_FORMULA 4 //布尔公式
double function1(double);
double function2(double);
double function3(double);
void Integral(int, double f(double), double, double); //矩形, 梯形逼近求定积分公式
void Trapezoid_Formula(double f(double), double a, double b); //递推梯形公式
void Simpson_Formula(double f(double), double a, double b); //递推辛普森公式
void Bool_Formula(double f(double), double a, double b); //递推布尔公式
void Formula(int, double f(double), double, double);
double function1(double x)
{
double y;
y = cos(x);
return y;
}
double function2(double x)
{
double y;
y=1/x;
// y = 2+sin(2 * sqrt(x));
return y;
}
double function3(double x)
{
double y;
y = exp(x);
return y;
}
int main()
{
printf("y=cos(x) , x=[0, 1]\n");
printf("Rectangle : "); Integral(RECTANGLE, function1, 0, 1);
printf("Trapezoid : "); Integral(TRAPEZOID, function1, 0, 1);
Formula(TRAPEZOID_FORMULA, function1, 0, 1);
Formula(SIMPSON_FORMULA, function1, 0, 1);
Formula(BOOL_FORMULA, function1, 0, 1);
printf("\n");
printf("y=1/x , x=[1, 5]\n");
printf("Rectangle : "); Integral(RECTANGLE, function2, 1, 5);
printf("Trapezoid : "); Integral(TRAPEZOID, function2, 1, 5);
Formula(TRAPEZOID_FORMULA, function2, 1, 5);
Formula(SIMPSON_FORMULA, function2, 1, 5);
Formula(BOOL_FORMULA, function2, 1, 5);
printf("\n");
printf("y=exp(x) , x=[0, 1]\n");
printf("Rectangle : "); Integral(RECTANGLE, function3, 0, 1);
printf("Trapezoid : "); Integral(TRAPEZOID, function3, 0, 1);
Formula(TRAPEZOID_FORMULA, function3, 0, 1);
Formula(SIMPSON_FORMULA, function3, 0, 1);
Formula(BOOL_FORMULA, function3, 0, 1);
return 0;
}
//积分:分割-近似-求和-取极限
//利用黎曼和求解定积分
void Integral(int method, double f(double),double a,double b) //f(double),f(x), double a,float b,区间两点
{
int i;
double x;
double sum = 0; //矩形总面积
double h; //矩形宽度
double t = 0; //单个矩形面积
h = (b-a) / DIV;
for(i=1; i <= DIV; i++)
{
x = a + h * i;
switch(method)
{
case RECTANGLE :
t = f(x) * h;
break;
case TRAPEZOID :
t = (f(x) + f(x - h)) * h / 2;
break;
default:
break;
}
sum = sum + t; //各个矩形面积相加
}
printf("the value of function is %.10f\n", sum);
}
//递推梯形公式
void Trapezoid_Formula(double f(double), double a, double b)//递推梯形公式
{
unsigned int i, j, M, J=1;
double h;
double F_2k_1 = 0;
double T[32];
double res = 0;
T[0] = (f(a) + f(b)) * (b-a)/2;
for(i=0; i<ITERATE; i++)
{
F_2k_1 = 0;
J = 1;
for(j=0; j<=i; j++) //区间划分
J <<= 1; //2^J
h = (b - a) /J; //步长
//M = J/2; //2M表示将积分区域划分的块数
for(j=1; j <= J; j += 2) //
{
F_2k_1 += f(a + j*h); //f(2k-1)累加和
}
T[i+1] = T[i]/2 + h * F_2k_1; //递推公式
res = T[i+1];
if(fabs(T[i+1] - T[i]) < EPS)
if(i < 3) continue;
else break;
}
printf("Trapezoid_Formula : the value of function is %.10f\n", res);
//return res;
}
//辛普森公式
void Simpson_Formula(double f(double), double a, double b) //辛普森公式
{
unsigned int i, j, M, J=1;
double h;
double F_2k_1 = 0;
double T[32], S[32];
double res_T = 0, res_S = 0;
T[0] = (f(a) + f(b)) * (b-a)/2;
for(i=0; i<ITERATE; i++)
{
F_2k_1 = 0;
J = 1;
for(j=0; j<=i; j++) //区间划分
J <<= 1; //2^J
h = (b - a) /J; //步长
//M = J/2; //2M表示将积分区域划分的块数
for(j=1; j <= J; j += 2) //
{
F_2k_1 += f(a + j*h); //f(2k-1)累加和
}
T[i+1] = T[i]/2 + h * F_2k_1; //递推梯形公式
res_T = T[i+1];
if(fabs(T[i+1] - T[i]) < EPS)
if(i < 3) continue;
else break;
}
for(i=1; i<ITERATE; i++)
{
S[i] = (4 * T[i] - T[i-1]) / 3; //递推辛普森公式
res_S = S[i];
if(fabs(S[i] - S[i-1]) < EPS)
if(i < 3) continue;
else break;
}
printf("Simpson_Formula : the value of function is %.10f\n", res_S);
//return res_S;
}
//布尔公式
void Bool_Formula(double f(double), double a, double b) //布尔公式
{
unsigned int i, j, M, J=1;
double h;
double F_2k_1 = 0;
double T[32], S[32], B[32];
double res_T = 0, res_S = 0, res_B;
T[0] = (f(a) + f(b)) * (b-a)/2;
for(i=0; i<ITERATE; i++)
{
F_2k_1 = 0;
J = 1;
for(j=0; j<=i; j++) //区间划分
J <<= 1; //2^J
h = (b - a) /J; //步长
//M = J/2; //2M表示将积分区域划分的块数
for(j=1; j <= J; j += 2) //
{
F_2k_1 += f(a + j*h); //f(2k-1)累加和
}
T[i+1] = T[i]/2 + h * F_2k_1; //递推公式
res_T = T[i+1];
if(fabs(T[i+1] - T[i]) < EPS)
if(i < 3) continue;
else break;
}
for(i=1; i<ITERATE; i++)
{
S[i] = (4 * T[i] - T[i-1]) / 3; //递推辛普森公式
res_S = S[i];
if(fabs(S[i] - S[i-1]) < EPS)
if(i < 3) continue;
else break;
}
for(i=1; i<ITERATE; i++)
{
B[i] = (16 * S[i] - S[i-1]) / 15; //递推布尔公式
res_B = B[i];
if(fabs(B[i] - B[i-1]) < EPS)
if(i < 3) continue;
else break;
}
printf("Bool_Formula : the value of function is %.10f\n", res_B);
//return res_B;
}
//采用二分区间的方法迭代,实际运行速度很慢
void Formula(int method, double f(double), double a, double b)//递推梯形公式
{
unsigned int i, j, M, J=1;
double h;
double F_2k_1 = 0;
double T[32], S[32], B[32];
double res_B = 0, res_S = 0, res_T = 0;
T[0] = (f(a) + f(b)) * (b-a)/2;
for(i=0; i<ITERATE; i++)
{
F_2k_1 = 0;
J = 1;
for(j=0; j<=i; j++) //区间划分
J <<= 1; //2^J
h = (b - a) /J; //步长
//M = J/2; //2M表示将积分区域划分的块数
for(j=1; j <= J; j += 2) //
{
F_2k_1 += f(a + j*h); //f(2k-1)累加和
}
T[i+1] = T[i]/2 + h * F_2k_1; //递推公式
res_T = T[i+1];
if(fabs(T[i+1] - T[i]) < EPS)
if(i < 3) continue;
else break;
}
switch(method)
{
default :
case TRAPEZOID_FORMULA :
printf("Trapezoid_Formula : the value of function is %.10f\n", res_T);
//return res_T;
break;
case SIMPSON_FORMULA :
for(i=1; i<ITERATE; i++)
{
S[i] = (4 * T[i] - T[i-1]) / 3;
res_S = S[i];
if(fabs(S[i] - S[i-1]) < EPS)
if(i < 3) continue;
else break;
}
printf("Simpson_Formula : the value of function is %.10f\n", res_S);
//return res_S;
break;
case BOOL_FORMULA :
for(i=1; i<ITERATE; i++)
{
S[i] = (4 * T[i] - T[i-1]) / 3;
res_S = S[i];
if(fabs(S[i] - S[i-1]) < EPS)
if(i < 3) continue;
else break;
}
for(i=1; i<ITERATE; i++)
{
B[i] = (16 * S[i] - S[i-1]) / 15;
res_B = B[i];
if(fabs(B[i] - B[i-1]) < EPS)
if(i < 3) continue;
else break;
}
printf("Bool_Formula : the value of function is %.10f\n", res_B);
//return res_B;
break;
}
}