数值分析第四章实验(数值积分)

suxss / 2023-08-01 / 原文

第一题

实现龙贝格积分算法，取 \(n=8,16,32,64,128\) 利用梯形公式，辛普森公式，龙贝格公式计算该积分，并与精确解进行比较

\[ \int_{-0.5}^{0.5} \frac{\sin x}{x} \, dx \]

import numpy as np
import pandas as pd
from scipy.integrate import quadrature

龙贝格公式

def romberg(x_min, x_max, f, tol=1e-6):  # 龙贝格公式
    h = (x_max - x_min)
    x = np.linspace(x_min, x_max, 2)
    y = f(x)
    data = [[(f(x_min) + f(x_max)) / 2], []]
    data[0].append(data[0][0] / 2 + h * f((x_min + x_max) / 2) / 2)
    data[1].append((4 * data[0][1] - data[0][0]) / 3)
    h /= 2
    x = np.linspace(x_min, x_max, 3)
    i = 2
    while (np.abs(data[i - 1][0] - data[i - 2][0]) > tol):
        pow_2_i = 4
        x_new = (x[:-1] + x[1:]) / 2
        y_new = f(x_new)
        h /= 2
        x = np.linspace(x_min, x_max, pow_2_i + 1)
        data[0].append(data[0][i - 1] / 2 + h * np.sum(y_new))
        data.append([])
        pow_4_j = 4
        for j in range(i):
            data[j + 1].append(
                pow_4_j * data[j][i - j] / (pow_4_j - 1) - data[j][i - j - 1] / (pow_4_j - 1))
            pow_4_j <<= 2
        i += 1
        pow_2_i << 1
    return data[-1][0], data

梯形公式

def trapezium(x_min, x_max, f, n):  # 梯形公式
    x = np.linspace(x_min, x_max, n + 1)
    y = f(x)
    return (x_max - x_min) * (np.sum(y) - (y[0] + y[-1]) / 2) / n

辛普森公式

def simpson(x_min, x_max, f, n):  # 辛普森公式
    w = np.array([4 if i % 2 else 2 for i in range(n + 1)])
    w[0] = w[-1] = 1
    x = np.linspace(x_min, x_max, n + 1)
    y = f(x)
    return np.sum(w * y) * (x_max - x_min) / (3 * n)

计算

def fun(x):
    if type(x) == np.ndarray :
        return np.array([fun(t) for t in x])
    else:
        return np.sin(x) / x if x != 0 else 1
		
n = [8, 16, 32, 64, 128]
I, _ = romberg(-0.5, 0.5, fun)
I_gauss = quadrature(fun, -0.5, 0.5)[0]
trapezium_values = [trapezium(-0.5, 0.5, fun, t) for t in n]
simpson_values = [simpson(-0.5, 0.5, fun, t) for t in n]
df = pd.DataFrame(data={"梯形公式": trapezium_values, "辛普森公式": simpson_values, "龙贝格公式": [I for t in n],
                        "自适应高斯": [I_gauss for t in n]})
df

\(n\)	梯形公式	辛普森公式	龙贝格公式	自适应高斯
\(8\)	\(0.985791\)	\(0.986215\)	\(0.986215\)	\(0.986215\)
\(16\)	\(0.986109\)	\(0.986215\)	\(0.986215\)	\(0.986215\)
\(32\)	\(0.986188\)	\(0.986215\)	\(0.986215\)	\(0.986215\)
\(64\)	\(0.986208\)	\(0.986215\)	\(0.986215\)	\(0.986215\)
\(128\)	\(0.986213\)	\(0.986215\)	\(0.986215\)	\(0.986215\)