<jupyter_start><jupyter_text>Дисциплина "Вычислительный практикум" Задание №5 Приближённое вычисление интеграла по квадратурным формулам наивысшей алгебраической степени точности Ковальчуков Александр, 223 группа Вариант №4  Постановка задачиТребуется приближенно вычислить интеграл$$\int_a^b \varphi(x)dx$$где подынтегральная функция может не быть достаточно гладкой на промежутке интегрирования,но представима в виде $\varphi(x) = \rho(x) f(x)$, где $\rho(x)$ содержит особенности, а $f(x)$ являетсядостаточно гладкой на $[a, b]$Далее будем рассматривать квадратурную формулу вида:$$\int_a^b \rho(x) f(x) dx \approx \sum_{k=1}^n A_k f(x_k) \label{eq:f1} \tag{1}$$где $x_k$ - узлы квадратурной формулы, а $A_k$ - коэффициенты$\rho(x)$ - называется весовой функцией, и для неё должны существовать моментывесовой функции:$$|\mu_k| = |\int_a^b \rho(x) x^k dx| < \infty, \;\; k = 0, 1, 2, \dots n$$Квадратурная формула будет интерполяционной, если$$A_k = \int_a^b \rho(x) \frac{\omega(x)}{(x - x_k) \omega'(x_k)} dx$$  Теорема 1 Для того, чтобы квадратурная формула $\eqref{eq:f1}$ была точна для любого многочлена степенине выше $2n - 1$, необходимо и достаточно, чтобы - узлы  $x_1, x_2, \dots, x_n$ являлись корнями ортогонального относительно веса $\rho(x)$ и отрезка $[a, b]$ многочлена $\omega(x) = (x - x_1)(x - x_2)\dots(x - x_n)$ - формула $\eqref{eq:f1}$ была интерполяционной  Теорема 2 Пусть отрезок интегрирования $[a, b]$ конечен. Если функция $f(x)$ имеет непрерывную на $[a, b]$производную порядка $2n$, то существует точка $\nu \in [a, b]$, такая что погрешность квадратурнойформулы $\eqref{eq:f1}$ гауссова типа имеет представление$$R_n(f) = \frac{f^{2n}(\mu)}{(2n)!} \int_a^b \rho(x) \omega^2(x)dx$$Код программы написан ня языке `python` с использованием интерактивной среды `Jupyter notebook`.<jupyter_code>import math
import pandas as pd<jupyter_output><empty_output><jupyter_text>Задача: Различными методами вычислить интеграл:$$\int_0^1 \frac{e^x}{\sqrt[3]x}$$<jupyter_code>I = 2.3435910933259677<jupyter_output><empty_output><jupyter_text>Точное значение интеграла: $$\int_0^1 \frac{e^x}{\sqrt[3]x} = \sqrt[3]{-1}(\Gamma(\frac{2}{3}, -1) - \Gamma(\frac{2}{3})) \approx 2.3435910933259677$$  Формула средних прямоугольников с тремя узлами<jupyter_code>def middle(a, b, func, w, m):
    h = (b - a) / m
    alpha = a + h/2
    return h * sum(w(alpha + k * h) * func(alpha + k * h) for k in range(m))

f = lambda x: math.e**x / x**(1/3)
w = lambda x: 1
M = 3
A = 0
B = 1
I_middle = middle(A, B, f, w, M)
print(I_middle)
print(abs(I_middle - I))<jupyter_output>2.2230282386930704
0.1205628546328974<jupyter_text>Значение интеграла по формуле средних прямоугольников с тремя узлами: $$\int_0^1 \frac{e^x}{\sqrt[3]x} \approx 2.2230282386930704$$  Величина абсолютной фактической погрешности:  0.120561561606013  Интерополяционная формула с весом $\frac{1}{\sqrt[3]{x}}$ по узлам $\frac{1}{6}, \; \frac{1}{2}, \; \frac{5}{6}$<jupyter_code>from sympy import Symbol, integrate, lambdify, solve, re
x = Symbol('x')

def definite_integral(func, a, b):
    indefinite_integral = lambdify(x, integrate(func, x))
    return indefinite_integral(b) - indefinite_integral(a)

A = 0
B = 1
x1 = 1 / 6
x2 = 1 / 2
x3 = 5 / 6
w = 1 / x**(1/3)
f =  lambda x: math.e**x
omega = (x - x1) * (x - x2) * (x - x3)
d_omega = lambda x: 3 * x**2 - 2 * (x1 + x2 + x3) * x + (x1 * x2 + x2 * x3 + x1 * x3)

mu0, mu1, mu2 = map(definite_integral, [w * x**k for k in range(3)], [A]*3, [B]*3)

A1 = 1 / d_omega(x1) * (mu2 - (x2 + x3) * mu1 + x2 * x3 * mu0)
A2 = 1 / d_omega(x2) * (mu2 - (x1 + x3) * mu1 + x1 * x3 * mu0)
A3 = 1 / d_omega(x3) * (mu2 - (x1 + x2) * mu1 + x1 * x2 * mu0)

I_interpolate = A1 * f(x1) + A2 * f(x2) + A3 * f(x3)
print(I_interpolate)
print(abs(I - I_interpolate))
print(f'x1: {x1}, x2: {x2}, x3: {x3}')
print(f'A1: {A1}, A2: {A2}, A3: {A3}')<jupyter_output>2.3459717130858713
0.002380619759903446
x1: 0.16666666666666666, x2: 0.5, x3: 0.8333333333333334
A1: 0.9000000000000001, A2: 0.14999999999999966, A3: 0.44999999999999996<jupyter_text>Значение интеграла по интерополяционной формуле с заданными узлами: $$\int_0^1 \frac{e^x}{\sqrt[3]x} \approx 2.3459717130858713 $$  Величина абсолютной фактической погрешности:  0.00238191278678812  Формула Гаусса с двумя узлами<jupyter_code>A = 0
B = 1
w = lambda x: 1
f = lambda x: math.e**x / x**(1/3)
t = 1 / 3**(1/2)

I_Gauss = (B - A) / 2 * (f((B - A) / 2 * - t + (B + A) / 2) + f((B - A) / 2 * t + (B + A) / 2))
print(I_Gauss)
print(abs(I - I_Gauss))
print("x1", (B - A) / 2 * - t + (B + A) / 2)
print("A1", (B - A) / 2)
print("x2", (B - A) / 2 * t + (B + A) / 2)
print("A2", (B - A) / 2)<jupyter_output>2.2278071408202447
0.11578395250572315
x1 0.21132486540518708
A1 0.5
x2 0.7886751345948129
A2 0.5<jupyter_text>Значение интеграла по формуле Гаусса: $$\int_0^1 \frac{e^x}{\sqrt[3]x} \approx 2.2278071408202447 $$  Величина абсолютной фактической погрешности:  0.11578395250572315  Формула типа Гаусса с двумя узлами<jupyter_code>import numpy as np
from sympy import Symbol, integrate, lambdify, solve, re
x = Symbol('x')

def definite_integral(func, a, b):
    indefinite_integral = lambdify(x, integrate(func, x))
    return indefinite_integral(b) - indefinite_integral(a)

def solve_kramer_rule(A, b):
    n, m = A.shape
    if n == m == len(b):
        solution = [0] * n
        det_A = np.linalg.det(A)
        for i in range(n):
            B = A.copy()
            B[:, i] = b
            solution[i] = np.linalg.det(B) / det_A
        return solution
    else:
        raise ValueError("matrix A and list b must have the same length")


A = 0
B = 1
w = 1 / x**(1/3)
f =  lambda x: math.e**x

mu0, mu1, mu2, mu3 = map(definite_integral, [w * x**k for k in range(4)], [A]*4, [B]*4)


A = np.array([[mu1, mu0], [mu2, mu1]])
b = np.array([-mu2, -mu3])

[a1, a2] = solve_kramer_rule(A, b)

print(a1, a2)

array = list(map(re, solve(x**2 + a1 * x + a2)))
# print(array)
[x1, x2] = array
print(f'x1: {x1}, x2: {x2}')
omega = (x - x1) * (x - x2)
d_omega = lambda x: 2 * x - (x1 + x2)



A1 = 1 / d_omega(x1) * (mu1 - x2 * mu0)
A2 = 1 / d_omega(x2) * (mu1 - x1 * mu0)

print(f"Коэффициенты квадратурной формулы: A1 = {A1}, A2 = {A2}")
print("Нулевой момент:", mu0)
print("Сумма коэффициентов:",A1 + A2)

I_Gauss_type = A1 * f(x1) + A2 * f(x2)
print(I_Gauss_type)
print(abs(I - I_Gauss_type))<jupyter_output>-0.9090909090909111 0.11363636363636441
x1: 0.149627093977301, x2: 0.759463815113608
Коэффициенты квадратурной формулы: A1 = 0.884164078649987, A2 = 0.615835921350013
Нулевой момент: 1.5
Сумма коэффициентов: 1.50000000000000
2.34299053839243
0.000600554933535502
{'fraction_non_alphanumeric': 0.1682764865680652, 'fraction_numerical': 0.1154542710534259, 'mean_word_length': 4.121329211746523, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 6, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '606611', 'license': 'MIT', 'n_tokens_mistral': 3544, 'n_tokens_neox': 3363, 'n_words': 734, 'path': 'lab_5.ipynb', 'repo_name': 'Qvery-mm/NumericMethods', 'size': '6626'}