Как посчитать интеграл методом трапеции и Симпсона на Python

Question

Вычислить приближённо интеграл точностью до 0.01

по формуле левых и правых прямоугольников, найдя число интервалов разбиения;
по формуле трапеций, найдя число интервалов разбиения;
по формуле Симпсона, найдя число интервалов разбиения;
по формуле трапеций, с помощью метода двойного пересчета;
по формуле Симпсона, с помощью метода двойного пересчета.

Методом прямоугольников посчитал, а вот с трапецией и Симпсоном не получается

from math import *
import numpy as np


def work(f, a, b, n):
   print("\nТекущее число разбиений: ", n)
   h = (b - a) / float(n)
   print("Текущий шаг:", h)
   total = sum([f((a + (k * h))) for k in range(0, n)])
   result = h * total
   print("Текущий результат: ", result)
   return result


# Туть подъинтегральная функция
def f(x):
    return (pow(x, 2)) * (pow(sin(x), 2))


print("Используем формулу левых прямоугольников")
print("Интегрируемая функция: f(x) = x^2*sin x^2")
print("Точность: 0.01")

n = 2
a1 = work(f, 1, 1.6, n)
n *= 2
a2 = work(f, 1, 1.6, n)

while abs(a1 - a2) > 0.01:
    n *= 2
    a1 = work(f, 1, 1.6, n)
    n *= 2
    a2 = work(f, 1, 1.6, n)

print("\nОтвет:", a2, "\nКоличество разбиений:", n)


def simpson(left, right, n, function):
   h = (right - left) / (2 * n)

   tmp_sum = function(left) + function(right)

   for step in range(1, 2 * n):
       if step % 2 != 0:
           tmp_sum += 4 * function(left + step * h)
       else:
           tmp_sum += 2 * function(left + step * h)

   return tmp_sum * h / 3


def trapezian(left, right, n, function):
    h = (right - left) / (n)

    return (function(left) + function(right) +
        sum(function(left + step * h) for step in range(1, n))) * h


print(simpson(0, 5, 15, f))
print(trapezian(0, 5, 30, f))

Answer 1

В формуле Симпсона умножается в конце на h / 6, а не на h / 3

Answer 2

Я не вижу понятно сформулированного условия. Вот пример вычисления методом трапеций для функции f(x) = x^2*sin x^2" на интервале 0:1 с точностью 0.01. Если это устроит буду разбираться с методом Симпсона.

import os

import numpy as np

def trapezoidal_rule(f, a, b, tol):
    n = 1
    integral_old = 0
    while True:
        x = np.linspace(a, b, n+1)
        y = f(x)
        h = (b - a) / n
        integral_new = (h / 2) * (y[0] + 2 * np.sum(y[1:-1]) + y[-1])
        if np.abs(integral_new - integral_old) < tol:
            break
        integral_old = integral_new
        n *= 2
    return integral_new

# Пример использования

print("-" * 50 + "\nПриближённое значение интеграла методом трапеций:\n" + "-" * 50)

f = lambda x: x**2 * np.sin(x**2)  # Функция, которую интегрируем
a = 0  # Нижний предел интегрирования
b = 1  # Верхний предел интегрирования
tol = 0.01  # Точность

result = trapezoidal_rule(f, a, b, tol)
print(f"Приближенное значение интеграла: {result}")

os.system("pause")

Теперь вычисление этого же методом Симпсона:

import os

import numpy as np

def simpson(f, a, b, tol):
    n = 2  # Начальное количество интервалов
    integral_old = 0  # Предыдущее значение интеграла
    while True:
        h = (b - a) / (2 * n)  # Шаг разбиения
        x = np.linspace(a, b, 2 * n + 1)  # Точки разбиения
        y = f(x)  # Значения функции в точках разбиения
        # Вычисляем новое значение интеграла по формуле Симпсона
        integral_new = h / 3 * (y[0] + 2 * np.sum(y[2:2*n:2]) + 4 * np.sum(y[1:2*n:2]) + y[2*n])
        # Проверяем условие остановки
        if np.abs(integral_new - integral_old) < tol:
            break
        integral_old = integral_new  # Обновляем предыдущее значение интеграла
        n *= 2  # Увеличиваем количество интервалов
    return integral_new

# Пример использования

print("-" * 50 + "\nПриближённое значение интеграла методом Симпсона:\n" + "-" * 50)

f = lambda x: x**2 * np.sin(x**2)
a = 0  # Нижний предел интегрирования
b = 1  # Верхний предел интегрирования
tol = 0.01  # Заданная точность
result = simpson(f, a, b, tol)  # Вычисляем интеграл
print(f"Приближенное значение интеграла: {result}")  # Выводим результат

os.system("pause")

БЛОГ НА HUSL

Как посчитать интеграл методом трапеции и Симпсона на Python

Ответы (2 шт):