基于复化辛卜生求积公式的变步长求积算法

1.辛卜生求积公式

对于I∗=∫abf(x)dxI^*=\int_{a}^{b}{f(x)dx}I∗=∫abf(x)dx ,差值型求积公式就是构造多项式Pn(x)P_n(x)Pn(x),使I∗≈I=I∗=∫abf(x)dxI^*\approx I= I^*=\int_{a}^{b}{f(x)dx}I∗≈I=I∗=∫abf(x)dx。
如果以a,c=a+b2,ba,c=\frac{a+b}{2},ba,c=2a+b,b为三个差值节点，构造二次插值多项式
P2(X)=(x−c)(x−b)(a−c)(a−b)f(a)+(x−a)(x−b)(c−a)(c−b)f(c)+(x−a)(x−c)(b−a)(b−c))f(b)P_2(X)=\frac{(x-c)(x-b)}{(a-c)(a-b)}f(a)+\frac{(x-a)(x-b)}{(c-a)(c-b)}f(c)+\frac{(x-a)(x-c)}{(b-a)(b-c))}f(b)P2(X)=(a−c)(a−b)(x−c)(x−b)f(a)+(c−a)(c−b)(x−a)(x−b)f(c)+(b−a)(b−c))(x−a)(x−c)f(b)
则可以推出I=∫abP2(x)dx==λ0f(a)+λ1f(c)+λ2f(b)I=\int_a^bP_2(x)dx==\lambda_0f(a)+\lambda_1f(c)+\lambda_2f(b)I=∫abP2(x)dx==λ0f(a)+λ1f(c)+λ2f(b)
其中λ0=∫ab(x−c)(x−b)(a−c)(a−b)dx=16(b−a)\lambda_0=\int_a^b\frac{(x-c)(x-b)}{(a-c)(a-b)}dx=\frac{1}{6}(b-a)λ0=∫ab(a−c)(a−b)(x−c)(x−b)dx=61(b−a)λ1=∫ab(x−a)(x−b)(c−a)(c−b)dx=46(b−a)\lambda_1=\int_a^b\frac{(x-a)(x-b)}{(c-a)(c-b)}dx=\frac{4}{6}(b-a)λ1=∫ab(c−a)(c−b)(x−a)(x−b)dx=64(b−a)λ2=∫ab(x−a)(x−c)(b−a)(b−c)dx=16(b−a)\lambda_2=\int_a^b\frac{(x-a)(x-c)}{(b-a)(b-c)}dx=\frac{1}{6}(b-a)λ2=∫ab(b−a)(b−c)(x−a)(x−c)dx=61(b−a)
称为辛卜生求积公式。
I∗=∫abf(x)dx≈b−a6[f(a)+4f(c)+f(b)]I^*=\int_{a}^{b}{f(x)dx}\approx \frac{b-a}{6}[f(a)+4f(c)+f(b)]I∗=∫abf(x)dx≈6b−a[f(a)+4f(c)+f(b)]
如果把[a,b]区间分为步长为h=b−anh=\frac{b-a}{n}h=nb−a的子区间[xk,xk+1][x_k , x_{k+1}][xk,xk+1]，再将子区间的定积分加起来得到整个区间的定积分近似值，该公式则称为复化辛卜生公式。

取步长为h=b−anh=\frac{b-a}{n}h=nb−a,将[a,b]分为n等分，分割点为xk=a+kh,k=0,1,2,3⋯n;x_k=a+kh,k=0,1,2,3\cdots n;xk=a+kh,k=0,1,2,3⋯n;
由辛卜生公式得：
∫abf(x)dx≈∑k=0n−1h6[f(xk)+4f(xk+xk+12)+f(xk+1)]\int_a^bf(x)dx\approx\sum_{k=0}^{n-1} \frac{h}{6} \left[ f(x_k)+4f(\frac{x_k+x_{k+1}}{2})+f(x_{k+1}) \right]∫abf(x)dx≈k=0∑n−16h[f(xk)+4f(2xk+xk+1)+f(xk+1)]
=h6[f(a)+4∑k=0n−14f(xk+xk+12)+2∑k=1n−1f(xk)+f(b)]=\frac{h}{6}\left[ f(a)+4\sum_{k=0}^{n-1}4f(\frac{x_k+x_{k+1}}{2})+2\sum_{k=1}^{n-1}f(x_k)+f(b)\right]=6h[f(a)+4k=0∑n−14f(2xk+xk+1)+2k=1∑n−1f(xk)+f(b)]
c语言程序实现如下：
例子∫12x2dx\int_1^2 x^2dx∫12x2dx：

#include <iostream>
double f(double x){return x*x;
}
double simpson(double (*f)(double x), double a, double b,int n){double result=0;double step=(b-a)/n;for (int i = 0; i <n ; ++i) {result+=f(a+step*i)+4*f(a+step*(i+1/2))+f(a+step*(i+1));}result=result*(step/6);return result;
}
int main(){double x;x=simpson(f,1,2,100000);printf("%lf",x);return 0;
}

结果：
2.333323
Process finished with exit code 0