数值分析作业-数值积分

HW

1. 确定数值积分 ${\int }_{-1}^{1}f(x)dx\approx \frac{2}{3}[f(-1)+f(0)+f(1)]$ 的迭代精度
2. 求 $A,x_0,x_1,x_2$ ，使得 ${\int }_{-1}^{1}f\left(x\right)dx\approx A[f\left(x_0\right)+f\left(x_1\right)+f\left(x_2\right)]$ 的代数精度至少为 3
3. 使用 2 点 Gauss 积分计算 ${\int }_0^1{x}^{3}dx$
4. 编程: 使用复合中心法则、复合梯形法则、复合 Simpson 法则、复合 3 点 Gauss 积分近似定积分 ${\int }_0^{1}x{e}^{x}dx$ ，与正确积分比较并报告 1、2、4、8、16、32 个子区间的误差。
5. 编程：使用 Romberg 积分近似定积分 ${\int }_0^{1}x{e}^{x}dx$ ,计算 $R_{jk}(1⩽j⩽5,1⩽k⩽j)$ 的误差.

1.

令 $f(x)=x$ , 积分计算得到上式左右两端成立，所以迭代精度至少为一。 令 $f(x)=x^2$ 积分计算得到两端等号不成立，所以精度为 1

2.

构造求解矩阵：

$$\left(\begin{array}{ccc}1& 1& 1\\ x_0& x_1& x_2\\ x_0^2& x_1^2& x_2^2\\ x_0^3& x_1^3& x_2^3\end{array}\right)\left(\begin{array}{c}A\\ A\\ A\end{array}\right)=\left(\begin{array}{c}2\\ 0\\ \frac{2}{3}\\ 0\end{array}\right)$$

解得： $A=\frac{2}{3},x_0=-\frac{1}{\sqrt{2}},x_1=0,x_2=\frac{1}{\sqrt{2}}$

3.

$$\begin{array}{rl}{\int }_0^1{x}^{3}dx& =\frac{1}{2}{\int }_{-1}^{1}f(\frac{t+1}{2})dt\\ & =\frac{1}{2}[(\frac{1-\sqrt{1/3}}{2}{)}^3+(\frac{1+\sqrt{1/3}}{2}{)}^3]\\ =0.250000& \end{array}$$

4.

程序运行截图

5

运行结果

| 1.3591409142295225 | 1.0917507747897934 | 1.0026207283098838 | 1.0230644790527572 | 1.0001690471404119 | 1.0000056017291137 | 1.0057741073678197 | 1.0000106501395072 | 1.0000000903394468 | 1.0000000028570712 | 1.0014440270677076 | 1.0000006669676702 | 1.0000000014228811 | 1.000000000011507 | 1.000000000000348|

1	2	3	4	5
1.3591409142295225
1.0917507747897934	1.0026207283098838
1.0230644790527572	1.0001690471404119	1.0000056017291137
1.0057741073678197	1.0000106501395072	1.0000000903394468	1.0000000028570712
1.0014440270677076	1.0000006669676702	1.0000000014228811	1.000000000011507	1.000000000000348

Appendix

3. Code

import math
 
# 编程: 使用复合中心法则、复合梯形法则、复合 Simpson 法则、复合 3 点 Gauss 积分近似定积分 $\int_{0}^{1} x e^{x} dx$ ，与正确积分比较并报告 1、2、4、8、16、32 个子区间的误差。
 
 
def Composite_Center(a, b, n):
    h = (b - a) / n
    x = [a + i * h for i in range(n + 1)]
    y = [math.exp((x[i] + x[i + 1]) / 2) for i in range(n)]
    return h * sum(y[i] for i in range(1, n))
 
 
def Composite_Trapezoid(a, b, n):
    h = (b - a) / n
    x = [a + i * h for i in range(n + 1)]
    y = [math.exp(x[i]) for i in range(n + 1)]
    return h / 2 * (y[0] + 2 * sum(y[i] for i in range(1, n)) + y[n])
 
 
def Composite_Simpson(a, b, n):
    h = (b - a) / n
    x = [a + i * h for i in range(n + 1)]
    y = [math.exp(x[i]) for i in range(n + 1)]
    yy = [math.exp((x[i] + x[i + 1]) / 2) for i in range(n)]
    return (
        h
        / 6
        * (
            y[0]
            + y[n]
            + 2 * sum(y[i] for i in range(1, n))
            + 4 * sum(yy[i] for i in range(n))
        )
    )  # noqa: E501
 
 
def gauss_legendre_3points(f, a, b):
    # Gauss-Legendre 3-points quadrature formula
    # weights and nodes
    w = [5 / 9, 8 / 9, 5 / 9]
    x0 = [-math.sqrt(3 / 5), 0, math.sqrt(3 / 5)]
    # change of variable
    tran = lambda x: 0.5 * (b - a) * x + 0.5 * (b + a)  # noqa: E731
    # transformed function
    g = lambda x: f(tran(x))  # noqa: E731
    # Gaussian quadrature formula
    Q = sum([w[i] * g(x0[i]) for i in range(3)])
    # change of variable in the integral
    op = 0.5 * (b - a) * Q
    return op
 
 
def composite_gauss_legendre_3points(f, a, b, n):
    # Composite Gauss-Legendre 3-points quadrature formula
    dx = (b - a) / n
    x = [a + i * dx for i in range(n + 1)]
    # subintervals
    return sum([gauss_legendre_3points(f, x[i], x[i + 1]) for i in range(n)])
 
 
def main():
    f = lambda x: x * math.e**x  # function to integrate  # noqa: E731
    a, b = 0, 1  # limits of integration
    n = [1, 2, 4, 8, 16, 32]  # number of subintervals
    for _ in n:
        op = composite_gauss_legendre_3points(f, a, b, _)
        cc = Composite_Center(a, b, _)
        cs = Composite_Simpson(a, b, _)
        ct = Composite_Trapezoid(a, b, _)
        print("________________________________________________________")
        print(f"Number of subintervals: {_}")
        print(f"Composite_Center:{(1-cc):.6E}")
        print(f"Composite_Trapezoid:{(1-ct):.6E}")
        print(f"Composite_Simpson:{(1-cs):.6E}")
        print(f"Composite_Gauss-Legendre:{(1-op):.6E}")
 
 
if __name__ == "__main__":
    main()
 

5.

import math
 
 
def Romberg(f, a, b, n):
    R = [[0 for _ in range(n + 1)] for _ in range(n + 1)]
    h = b - a
    R[0][0] = 0.5 * h * (f(a) + f(b))
    powerOf2 = 1
    for i in range(1, n + 1):
        h = 0.5 * h
        R[i][0] = 0.5 * R[i - 1][0] + h * sum(f(a + (2 * j - 1) * h) for j in range(1, powerOf2 + 1))
        powerOf2 *= 2
        for k in range(1, i + 1):
            R[i][k] = R[i][k - 1] + (R[i][k - 1] - R[i - 1][k - 1]) / (4 ** k - 1)
    return R
 
 
def main():
    f = lambda x: x * math.exp(x)
    a, b = 0, 1
    n = 5
    R = Romberg(f, a, b, n)
    print(R)
 
 
if __name__ == "__main__":
    main()