梯形求积公式 和 复合梯形求积公式 Matlab 实现



梯形求积公式


  
仅使用区间两点x1,f(x1),x2,f(x2)

  组成的梯形面积S代替∫x2x1f(x)dx的近似方法 

∫x2x1f(x)dx≈S=x2−x12×(f(x2)−f(x1))



复合梯形求积公式


  
将求积区间[a,b]分为n个区间，每个区间步长为h(h=b−an)在每个区间求梯形积分 


  Si为第i个梯形的面积 
∫x2x1f(x)dx≈∑i=0nSi=h2∫(f(a)+∑k=1n−1f(xk)+f(b))



说明


  
上述公式是我用mathjax写的，如有错误请联系我修正  
敬请指正  

  概述省略了部分推导过程，请查阅详细推导资料




Matlab 实现代码



梯形求积公式


  
将该函数存为m文件




function res = Trapezium(f,a,b)
    format long;
    if b < a
        c = b;
        b = a;
        a = c;
    end
    res = (b-a)*(f(a)+f(b))/2;


  
调用下面语句测试函数




f = inline('sin(x)','x')
Trapezium(f,0,pi/2)



复合梯形求积公式


  
将该函数存为m文件




function res = ComTrapezium(f,n,a,b)
    format long;
    if b < a
        c = b;
        b = a;
        a = c;
    end
    h = (b-a)/n;
    d = f(a);
    for i = a+h:h:b-h 
        d = d + (2 * f(i));
    end
    d = d + f(b);
    res = (d * h / 2);


  
调用下面语句测试函数




f = inline('sin(x)','x')
ComTrapezium(f,4,0,pi/2)

【C语言基础】利用复合梯形求积公式计算定积分
一、复合梯形求积公式
这是数值分析中一种求解定积分的近似方法。适用于被积函数的原函数不能用初等函数表示的情况。
基本思路
将被积函数 f（x）与x轴围成的区域分成n个梯形，把n个梯形面积求和得到积分的近似值。若精度不满足需要，则可以将每个区间再等分一次，得到2n+1个等分区间，然后再求和，直到精度满足需要。
公式
$\int\limits_a^b f(x)dx \approx \frac{h}{2}[f(a)+f(b)+2\sum_{i=1}^{n-1} f(x_i)]$

其中：
n 为等分区间数
h=（b-a）/2 (积分步长)
C语言实现算法
以下面积分为例：
$\int\limits_0^11+e^-xsin4xdx$
首先写出计算被积函数值的函数
double f(double x)
{
    return 1 + exp(-x)*sin(4x);
}

然后写出求积分近似解的表达式
double caculate(double a,double b)
{
    double y1,y2,h;
    int n=1,i;
    h = a-b;
    y1 = h*(f(a)+f(b))/2;
    do
    {
       y2 = y1;
       n = n*2;
       h = (a-b)/n;
       y1 = f(a)+f(b);
       for(i=1;i<n:i++)
             y1 = y1 + f(a + i*h);
        y1 = h*y1/2;
    }while(fabs(y2-y1)>1.0e-5));
    return y2;
}

其中frab（）是用来计算绝对值的函数

实验内容


计算圆周率π，要求误差小于10-8。
	
用复合梯形公式和复合辛普森求积公式计算。  
 

 

​


复合梯度代码
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function T_n=fht(a,b,n) 
h=(b-a)/n; 
for k=0:n     
    x(k+1)=a+k*h;     
    if x(k+1)==0          
        x(k+1)=10^(-10);     
    end
end
T_1=h/2*(fx1(x(1))+fx1(x(n+1))); 
for i=2:n      
    F(i)=h*fx1(x(i)); 
end
T_2=sum(F); 
T_n=T_1+T_2;

  复合辛普森求积代码
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function S_n=S_P_S(a,b,n) 
h=(b-a)/n;
for k=0:n     
    x(k+1)=a+k*h;      
    x_k(k+1)=x(k+1)+1/2*h;     
    if (x(k+1)==0)|(x_k(k+1)==0)         
        x(k+1)=10^(-10);         
        x_k(k+1)=10^(-10); 
    end
end
S_1=h/6*(fx1(x(1))+fx1(x(n+1)));
for i=2:n      
     F_1(i)=h/3*fx1(x(i)); 
end
for j=1:n      
     F_2(j)=2*h/3*fx1(x_k(j)); 
end
S_2=sum(F_1)+sum(F_2); 
S_n=S_1+S_2;
                                     龙贝格和变步长公式代码
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function asd
format long;
   a=0;
   b=1;
   m=2;
   t(1)=0.5*(b-a)*(f(a)+f(b));
   t(2)=0.5*t(1)+0.5*(b-a)*f((a+b)/2);
   s(1)=(4*t(2)-t(1))/3;
   j=2;
   while abs (t(j)-t(j-1))>(0.5e-8/4),
       h=(b-a)/m;
       k=0:(m-1);
       j=j+1;
       t(j)=0.5*t(j-1)+0.5*h*sum(f(a+(k+1/2)*h));
       s(j-1)=(4*t(j)-t(j-1))/3;
       c(j-2)=(16*s(j-1)-s(j-2))/15;
       if j>3,
           r(j-3)=(64*c(j-2)-c(j-3))/63;
       end
       m=m*2;
   end
   fsimpson=4*s(end)
   rbg=4*r(end)
   bbb=4*c(end)
   j
                                  方程代码
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function y=fx1(x) 
y=4/(1+x^2);

[点击并拖拽以移动]
​

 

>> T_20=fht(0,1,1048576)

T_20 =

    3.1416

>> vpa(T_20,20)
 
ans =
 
3.1415926535896492311
 
>> 




S_1=S_P_S(0,1,2)

S_1 =

    3.2200

>>  vpa(S_1,7) 
 
ans =
 
3.22
 
>>  S_2=S_P_S(0,1,4)

S_2 =

    3.1518

>> vpa(S_2,7)
 
ans =
 
3.151849
 
>> S_3=S_P_S(0,1,8)

S_3 =

    3.1429

>> vpa(S_3,7)
 
ans =
 
3.14289
 
>> S_4=S_P_S(0,1,16)

S_4 =

    3.1418


bbb =

   3.1415955553589793


bbb =

   3.14159244666653589793

bbb =

   3.141592653589793

 

1.求积公式余项
1.1 定义
$R[f]=\int_a^b\!\!\!f(x)\mathrm{d}x-\sum_{k=0}^nA_kf(x_k)=Kf^{(m+1)}(\eta ),(1)$

$其中K为不依赖f(x)的待定参数,\eta \in (a,b).$$这个结果表明当f(x)是次数小于等于m的多项式时,$

$由于f^{(m+1)}(x)=0,故此时R[f]=0,即求积公式\int_a^b\!f(x)\mathrm{d}x\approx\sum\limits_{k=0}^nA_kf(x_k)$

$精确成立.而当f(x)=x^{m+1}时,f^{(m+1)}(x)=(m+1)!,(1)式左端R_n(f)\neq 0,故可求得$

$\begin{aligned}
    K&=\frac{1}{(m+1)!}\left[\int_a^b\!\!\!x^{m+1}\mathrm{d}x-\sum_{k=0}^nA_kx^{m+1}\right]\\
    &=\frac{1}{(m+1)!}\left[\frac{1}{(m+2)}(b^{m+2}-a^{m+2})-\sum_{k=0}^nA_kx^{m+1}\right]
\end{aligned}$
1.2 Python实现求积公式余项
from sympy.abc import x
from sympy import integrate, exp, diff
import numpy as np
from math import factorial


def quadrature_reminder(a, b, equal_segment, m, intergrand, a_k: np.ndarray):
    """
    实现[a,b]代数精度为m的求积公式余项表达式
    :param a: 区间左端点即起点a
    :param b: 区间右端点即终点b
    :param equal_segment: 区间等分数n
    :param m: 代数精度m
    :param intergrand: 被积函数f(x)
    :param a_k: 权系数Ak数组
    :return: 求积公式余项表达式
    """
    step_h = (b - a) / equal_segment
    f_x = x ** (m + 1)
    f_x_k_array = np.array([f_x.subs(x, a + k * step_h) for k in range(equal_segment + 1)])
    return 1 / factorial(m + 1) * (1 / (m + 2) * (b ** (m + 2) - a ** (m + 2)) - np.sum(a_k * f_x_k_array)) * diff(
        intergrand, x, m + 1), a, b

2.牛顿-柯特斯公式
2.1 定义
$设将积分区间[a,b]分为n等分,步长h=\dfrac{b-a}{n},选取等距节点x_k=a+kh$$构造出的插值型求积公式$

$I_n=(b-a)\sum_{k=0}^nC_k^{(n)}f(x_k),$

$称为\textbf{牛顿-柯特斯}(\mathrm{Newtown-Cotes})\textbf{公式},式中C_k^{(n)}称为\textbf{柯特斯系数}.$

$根据求积系数A_k=\int_a^bl_k(x)\mathrm{d}x,k=0,1,\ldots,n.引进变换x=a+th,则有$

$C_k^{(n)}=\frac{h}{b-a}\int_0^n\prod_{j=0\\j\neq k}^n\frac{t-j}{k-j}\mathrm{d}t=\frac{(-1)^{n-k}}{nk!(n-k)!}\int_0^n\prod_{j=0\\j\neq k}^n(t-j)\mathrm{d}t.$

$特别地，当n=2时,相应的求积公式称为\textbf{辛普森}(\mathrm{Simpson})\textbf{公式}:$

$S=\frac{b-a}{6}\left[f(a)+4f(\frac{a+b}{2})+f(b)\right].$
2.2 Python实现牛顿-柯特斯公式
def newtown_cotes_integrate(equal_segment, interval_start, interval_end, symbol_t, f_x):
    """
    实现牛顿-柯特斯(Newton-Cotes)插值型积分
    :param symbol_t: 引进变换x=a+t*h后的积分变量
    :param equal_segment: 区间等分数n
    :param interval_start: 区间左端点即起点
    :param interval_end: 区间右端点即终点
    :param f_x: 被积函数 intergrand
    :return:牛顿-柯特斯(Newton-Cotes)插值积分
    """
    step_h = (interval_end - interval_start) / equal_segment
    f_k_array = np.array([f_x.subs(x, interval_start + k * step_h) for k in range(equal_segment + 1)])
    return (interval_end - interval_start) * np.sum(cotes_coefficient(equal_segment, symbol_t) * f_k_array)


def cotes_coefficient(equal_segment, symbol_t):
    """
    实现牛顿-柯特斯(Newton-Cotes)插值型求积公式的柯特斯系数
    :param equal_segment: 区间等分数n,其中1<=n<=7,n>7的牛顿-柯特斯公式计算不稳定，不使用.
    :param symbol_t: 引进变换x=a+t*h后的积分变量t
    :return: 柯特斯系数数组
    """
    if equal_segment not in range(1, 8):
        raise ValueError("Cotes coefficient must be an integer between 1 and 7")
    c_k_array = np.zeros(equal_segment + 1, dtype=np.float64)
    for k in range(equal_segment + 1):
        index = list(range(equal_segment + 1))
        index.pop(k)
        intergrand = np.prod(symbol_t - np.array(index))
        integral = integrate(intergrand, (symbol_t, 0, equal_segment))
        c_k_array[k] = (-1) ** (equal_segment - k) / (
                equal_segment * factorial(k) * factorial(equal_segment - k)) * integral
    return c_k_array

3.复合梯形公式
3.1 定义
$将区间[a,b]划分为n等份,分点x_k=a+kh,h=\dfrac{b-a}{n},k=0,1,\ldots,n,$$在每个子区间[x_k,x_{k+1}](k=0,1,\ldots,n-1)上$

$采用梯形公式\int_a^b\!f(x)\mathrm{d}x\approx\dfrac{b-a}{2}[f(a)+f(b)],则得$

$T_n=\frac{h}{2}\sum_{k=0}^{n-1}[f(x_k)+f(x_{k+1})]=\frac{h}{2}\left[f(a)+2\sum_{k=1}^{n-1}f(x_k)+f(b)\right],$

称为复合梯形公式
3.2 Python实现复合梯形公式
def com_trapz(a, b, equal_segment, intergrand):
    """
    实现复合梯形求积公式
    :param a: 区间左端点即起点a
    :param b: 区间右端点即终点b
    :param equal_segment: 区间等分数n
    :param intergrand: 被积函数f(x)
    :return: 复合梯形求积积分
    """
    step_h = (b - a) / equal_segment
    f_x_k_array = np.array([intergrand.subs(x, a + k * step_h) for k in range(1, equal_segment)])
    return step_h / 2 * (intergrand.subs(x, a) + 2 * np.sum(f_x_k_array) + intergrand.subs(x, b))

4.复合辛普森公式
4.1 定义
$将区间[a,b]划分为n等份,在每个子区间[x_k,x_{k+1}](k=0,1,\ldots,n-1)上$

$采用辛普森公式,若记x_{k+1/2}=x_k+\frac{1}{2}h,则得$

$\begin{aligned}
    S_n&=\frac{h}{6}\sum_{k=0}^{n-1}[f(x_k)+4f(x_{k+1/2})+f(x_{k+1})]\\
    &=\frac{h}{6}\left[f(a)+4\sum_{k=0}^{n-1}f(x_{k+1/2})+2\sum_{k=1}^{n-1}f(x_k)+f(b)\right],
\end{aligned}$

称为复合辛普森公式
4.2 Python实现复合辛普森公式
def com_simpson(a, b, equal_segment, intergrand):
    """
    实现复合辛普森求积公式
    :param a: 区间左端点即起点a
    :param b: 区间右端点即终点b
    :param equal_segment: 区间等分数n
    :param intergrand: 被积函数f(x)
    :return: 复合辛普森求积积分
    """
    step_h = (b - a) / equal_segment
    f_x_k1_array = np.array([intergrand.subs(x, a + (k+0.5) * step_h) for k in range(equal_segment)],
                            dtype=np.float64)
    f_x_k2_array = np.array([intergrand.subs(x, a + k * step_h) for k in range(1, equal_segment)],
                            dtype=np.float64)
    return step_h / 6 * (
            intergrand.subs(x, a) + 4 * np.sum(f_x_k1_array) + 2 * np.sum(f_x_k2_array) + intergrand.subs(x, b))

5.测试
if __name__ == '__main__':
    # 辛普森公式及其余项表达式测试成功，来源详见来源详见李庆扬数值分析第5版P135,e.g.4
    # 0.632333680003663
    inter_grand = exp(-x)
    print("辛普森公式积分为:{}".format(newtown_cotes_integrate(2, 0, 1, x, inter_grand)))
    print("辛普森公式为:{}".format(quadrature_reminder(a=0, b=1, equal_segment=2, m=3, intergrand=exp(-x), a_k=cotes_coefficient(2, x))[0]))
    # 复合梯形公式和复合辛普森公式测试成功，来源详见来源详见李庆扬数值分析第5版P135,e.g.2(1)
    f_x = x / (4 + x ** 2)
    print("复合梯形公式积分:{}".format(com_trapz(a=0, b=1, equal_segment=8, intergrand=f_x)))  # 0.111402354529548
    print("复合辛普森公式积分:{}".format(com_simpson(a=0, b=1, equal_segment=8, intergrand=f_x)))  # 0.111571813252631

6.运行结果