数值积分

一. 组合梯形公式（数值积分）

【问题描述】组合梯形公式求函数f(x)=2+sin(2*sqrt(x))的积分近似值。

【输入形式】在屏幕上依次输入积分上限、下限和等距子区间个数。

【输出形式】输出使用组合梯形公式求得的积分近似值。

【样例1输入】

1 6 10

【样例1输出】

8.19385457

【样例1说明】输入：积分上限a为1、下限b为6和等距子区间个数m为10。输出：积分近似值（保留小数点后8位有效数字）

【评分标准】根据输入得到的输出准确

二. 组合辛普森公式（数值积分）

【问题描述】组合辛普森公式求函数f(x)=2+sin(2*sqrt(x))的积分近似值。

【输入形式】在屏幕上依次输入积分上限、下限和等距子区间个数。

【输出形式】输出使用组合辛普森公式求得的积分近似值。

【样例1输入】

1 6 5

【样例1输出】

8.18301549

【样例1说明】输入：积分上限a为1、下限b为6和等距子区间个数m为5。输出：积分近似值（保留小数点后8位有效数字）

【评分标准】根据输入得到的输出准确

import math

def fun(x):
    return math.sin(2*math.sqrt(x))+2

# 复合梯度
def tx(a,b,m):
    h = (b-a)/m
    s = 0
    for i in range(1,m):
        x = a+h*i
        s = s+fun(x) #s=f1+f2+----+fm-1
    result = (h/2)*(fun(a)+s*2+fun(b))
    return result

#复合辛普森
def xps(a,b,m):
    h = (b-a)/(2*m)
    s1 = 0
    s2 = 0
    for i in range(1,m+1):
        x = a+ h*(2*i -1)
        s1 = s1+fun(x)
    s1 = 4*h*s1/3
    for i in range(1,m):
        x = a+ h*2*i
        s2 = s2+fun(x)
    s2 = 2 * h * s2 / 3
    result = (h/3)*(fun(a)+fun(b))+s1+s2
    return result
 """
 def xps(a,b,m):
    h = (b-a)/(2*m)
    s1 = 0
    for i in range(1,m+1):
        x = a+ h*2*i
        s1 = s1+fun(x)+4*fun(x-h)+fun(x-2*h)
    result = (h/3)*s1
    return result
"""


def main():
    a,b,m= input().split()
    a =int(a)
    b =int(b)
    m =int(m)
    print("%.8f"%xps(a, b, m))
if __name__ == '__main__':
    main()

三. 龙贝格积分（数值积分）

【问题描述】利用龙贝格积分公式计算函数f(x)=(x^2+x+1)cos(x)，在区间[0, pi/2]范围内的定积分近似值。

【输入形式】在屏幕上龙贝格积分表行的最大值。

【输出形式】龙贝格积分表。

【样例1输入】

6

【样例1输出】

[[0.78539816 0. 0. 0. ]

[1.72681266 2.04061749 0. 0. ]

[1.96053417 2.03844134 2.03829626 0. ]

[2.01879395 2.03821388 2.03819871 2.03819716]

[2.03334734 2.03819847 2.03819745 2.03819743]

[2.03698495 2.03819749 2.03819743 2.03819743]]

【样例1说明】输入：龙贝格积分表行的最大值n为6。输出：6行4列的龙贝格积分表。

【评分标准】根据输入得到的输出准确

import math
import numpy as np

def fun(x):
    return (x**2 + x + 1) * math.cos(x)
#龙贝格积分
def romber(a,b,n):
    # M为区间个数
    M = 1
    # h为区间长度
    h = b-a
    err = 1
    # j是 2**j=2M
    j = 0
    R = np.zeros([n, 4], dtype=np.double)
    R[0,0] = h*(fun(a)+fun(b))/2
    # 第一列
    while( j < n-1 ):
        h = h/2
        s = 0
        for k in range(1, M+1):
            x = a + h*(2*k-1)
            s = s + fun(x)
        R[j+1, 0] = R[j, 0]/2 + h*s
        j = j + 1
        M = 2*M
    #向下算
    for k in range(1,4):
        for j in range(k,n):
             R[j, k] = R[j, k-1] + (R[j, k-1]-R[j-1, k-1])/(4**k-1)
    print(R)

def main():
    n = int(input())
    a = 0
    b = math.pi/2
    romber(a, b, n)
if __name__ == '__main__':
    main()