Глава 3. МАТЕМАТИЧЕСКИЙ АНАЛИЗ (CALCULUS)
3.1. Определение производной через предел
\[ f'(x) = \lim_{h \to 0} \frac{f(x + h) - f(x)}{h} \]Объяснение: Производная функции определяется как предел отношения приращения функции к приращению аргумента, когда приращение аргумента \(h\) стремится к нулю. Она представляет мгновенную скорость изменения функции.
Пример: Для \(f(x) = x^2\), вычислить \(f'(x) = \lim_{h \to 0} \frac{(x+h)^2 - x^2}{h} = \lim_{h \to 0} \frac{x^2 + 2xh + h^2 - x^2}{h} = \lim_{h \to 0} \frac{2xh + h^2}{h} = \lim_{h \to 0} (2x + h) = 2x\).
Реализация (численное дифференцирование):
def derivative(f, x, h=1e-5):
"""Вычисляет численную производную функции f в точке x."""
return (f(x + h) - f(x)) / h
# Пример использования:
def my_func(x):
return x**2
deriv_at_3 = derivative(my_func, 3)
print(deriv_at_3) # Должно быть близко к 6
3.2. Степенное правило
\[ \frac{d}{dx} x^n = nx^{n-1} \]Объяснение: Степенное правило упрощает дифференцирование одночленов. Оно является основополагающим в матанализе и широко используется в вычислениях градиентов в МО.
Пример: Для \(f(x) = x^3\), \(f'(x) = 3x^{3-1} = 3x^2\).
Реализация:
def power_rule(n, x):
"""Применяет степенное правило дифференцирования."""
if n == 0:
return 0 # Производная константы
else:
return n * x**(n - 1)
# Пример
print(power_rule(3, 5)) # Выведет 3 * 5**2 = 75
3.3. Правило произведения
\[ \frac{d}{dx} [u(x)v(x)] = u'(x)v(x) + u(x)v'(x) \]Объяснение: Правило произведения вычисляет производную произведения двух функций. Оно критически важно для работы с мультипликативными зависимостями в МО.
Пример: Для \(f(x) = (x^2)(e^x)\), \(u(x)=x^2, v(x)=e^x\). \(u'(x)=2x, v'(x)=e^x\). Тогда \(f'(x) = (2x)(e^x) + (x^2)(e^x) = 2xe^x + x^2e^x\).
Реализация:
# Предполагается, что u, v, u_prime, v_prime - это функции
def product_rule(u, v, u_prime, v_prime, x):
"""Применяет правило произведения для дифференцирования."""
return u_prime(x) * v(x) + u(x) * v_prime(x)
# Пример использования:
def u(x): return x**2
def v(x): return np.exp(x)
def u_prime(x): return 2*x
def v_prime(x): return np.exp(x)
import numpy as np
deriv_product = product_rule(u, v, u_prime, v_prime, 1)
print(deriv_product) # 2*1*e^1 + 1^2*e^1 = 3e ~ 8.15
3.4. Правило частного
\[ \frac{d}{dx} \left[ \frac{u(x)}{v(x)} \right] = \frac{u'(x)v(x) - u(x)v'(x)}{[v(x)]^2} \]Объяснение: Правило частного вычисляет производную отношения двух функций. Оно необходимо для операций, включающих деление в моделях МО.
Пример: Для \(f(x) = \frac{x^2}{e^x}\), \(u(x)=x^2, v(x)=e^x\). \(u'(x)=2x, v'(x)=e^x\). Тогда \(f'(x) = \frac{(2x)(e^x) - (x^2)(e^x)}{(e^x)^2} = \frac{2xe^x - x^2e^x}{e^{2x}} = \frac{(2x - x^2)e^x}{e^{2x}} = \frac{2x - x^2}{e^x}\).
Реализация:
# Предполагается, что u, v, u_prime, v_prime - это функции
def quotient_rule(u, v, u_prime, v_prime, x):
"""Применяет правило частного для дифференцирования."""
v_x = v(x)
if v_x == 0:
raise ValueError("Деление на ноль в знаменателе")
return (u_prime(x) * v_x - u(x) * v_prime(x)) / (v_x**2)
# Пример использования (с функциями из предыдущего примера):
deriv_quotient = quotient_rule(u, v, u_prime, v_prime, 1)
print(deriv_quotient) # (2*1*e^1 - 1^2*e^1) / (e^1)^2 = e / e^2 = 1/e ~ 0.367
3.5. Цепное правило (Правило дифференцирования сложной функции)
\[ \frac{d}{dx} f(g(x)) = f'(g(x)) g'(x) \]Объяснение: Цепное правило вычисляет производную сложной функции. Оно широко используется в обратном распространении ошибки для обучения нейронных сетей.
Пример: Для \(f(x) = \sin(x^2)\). Пусть \(h(x) = \sin(x)\) и \(g(x) = x^2\). Тогда \(f(x) = h(g(x))\). \(h'(x) = \cos(x)\) и \(g'(x) = 2x\). По цепному правилу \(f'(x) = h'(g(x)) \cdot g'(x) = \cos(x^2) \cdot 2x\).
Реализация:
# f_prime - производная внешней функции
# g - внутренняя функция
# g_prime - производная внутренней функции
def chain_rule(f_prime, g, g_prime, x):
"""Применяет цепное правило."""
return f_prime(g(x)) * g_prime(x)
# Пример использования:
def g(x): return x**2
def g_prime(x): return 2*x
def f_prime(y): return np.cos(y) # Производная от sin(y)
import numpy as np
deriv_chain = chain_rule(f_prime, g, g_prime, np.pi/2)
# cos((pi/2)^2) * 2*(pi/2) = cos(pi^2/4) * pi
print(deriv_chain)
3.6. Производная логарифма
\[ \frac{d}{dx} \ln(x) = \frac{1}{x}, \quad x > 0 \]Объяснение: Производная функции натурального логарифма равна обратной величине ее аргумента. Часто используется в МО для оптимизации и логарифмических преобразований.
Пример: Для \(f(x) = \ln(x)\), \(f'(2) = \frac{1}{2}\).
Реализация:
import numpy as np
def log_derivative(x):
"""Вычисляет производную натурального логарифма."""
if x <= 0:
raise ValueError("Аргумент логарифма должен быть положительным")
return 1 / x
print(log_derivative(2)) # Выведет 0.5
3.7. Производная экспоненты
\[ \frac{d}{dx} e^x = e^x \]Объяснение: Экспоненциальная функция уникальна тем, что ее производная равна самой функции. Это свойство является ключевым в вычислениях градиентов и моделях экспоненциального роста в МО.
Пример: Для \(f(x) = e^x\), \(f'(2) = e^2\).
Реализация:
import numpy as np
def exp_derivative(x):
"""Вычисляет производную экспоненты e^x."""
return np.exp(x)
print(exp_derivative(2)) # Выведет e^2 ~ 7.389
3.8. Интеграл степенной функции
\[ \int x^n dx = \frac{x^{n+1}}{n + 1} + C, \quad n \neq -1 \]Объяснение: Интеграл степенной функции обобщает нахождение первообразной для одночленов. Это правило фундаментально в интегральном исчислении и применяется в МО для анализа функций стоимости.
Пример: Для \(f(x) = x^2\), \(\int x^2 dx = \frac{x^{2+1}}{2+1} + C = \frac{x^3}{3} + C\).
Реализация (нахождение первообразной):
import sympy
x_sym = sympy.symbols('x')
def power_integral_symbolic(n, x_var):
"""Находит символьное выражение первообразной степенной функции."""
if n == -1:
return sympy.log(x_var) # Первообразная 1/x это ln|x|, здесь для x>0
else:
return x_var**(n + 1) / (n + 1)
# Пример использования
integral_expr = power_integral_symbolic(2, x_sym)
print(integral_expr) # Выведет x**3/3
Примечание: Python код обычно вычисляет определенные интегралы численно, а не находит символьные первообразные. Пример выше использует SymPy для демонстрации правила.
3.9. Основная теорема матанализа (Формула Ньютона-Лейбница)
\[ \int_{a}^{b} f(x)dx = F(b) - F(a), \quad \text{где } F'(x) = f(x) \]Объяснение: Основная теорема матанализа связывает дифференцирование и интегрирование, утверждая, что определенный интеграл по интервалу равен разности значений первообразной на концах интервала.
Пример: Для \(f(x) = x^2\) на отрезке \([1, 3]\). Первообразная \(F(x) = \frac{x^3}{3}\). \(\int_{1}^{3} x^2 dx = F(3) - F(1) = \frac{3^3}{3} - \frac{1^3}{3} = \frac{27}{3} - \frac{1}{3} = \frac{26}{3}\).
Реализация (численное интегрирование):
from scipy.integrate import quad
def definite_integral(f, a, b):
"""Вычисляет определенный интеграл функции f от a до b."""
result, error_estimate = quad(f, a, b)
return result
# Пример использования:
def my_func(x):
return x**2
integral_value = definite_integral(my_func, 1, 3)
print(integral_value) # Должно быть близко к 26/3 ~ 8.667
3.10. Частные производные
\[ \frac{\partial f}{\partial x} = \lim_{h \to 0} \frac{f(x + h, y) - f(x, y)}{h}, \quad \frac{\partial f}{\partial y} = \lim_{h \to 0} \frac{f(x, y + h) - f(x, y)}{h} \]Объяснение: Частные производные измеряют скорость изменения функции многих переменных по одной переменной, при этом остальные переменные считаются константами. Они необходимы в оптимизации и градиентных методах МО.
Пример: Для \(f(x, y) = x^2 + y^2\), \(\frac{\partial f}{\partial x} = 2x\), \(\frac{\partial f}{\partial y} = 2y\).
Реализация (численное дифференцирование):
def partial_derivative(f, var_index, point, h=1e-5):
"""Вычисляет численную частную производную функции f."""
args1 = list(point)
args2 = list(point)
args1[var_index] += h
# Используем центральную разность для большей точности (опционально)
# args2[var_index] -= h
# return (f(*args1) - f(*args2)) / (2 * h)
# Прямая разность (как в определении)
return (f(*args1) - f(*point)) / h
# Пример использования:
def my_multi_func(x, y):
return x**2 + y**2
point = (1, 2)
partial_x = partial_derivative(my_multi_func, 0, point) # по x (индекс 0)
partial_y = partial_derivative(my_multi_func, 1, point) # по y (индекс 1)
print(f"df/dx at {point} ~ {partial_x}") # Должно быть близко к 2*1 = 2
print(f"df/dy at {point} ~ {partial_y}") # Должно быть близко к 2*2 = 4
3.11. Градиент
\[ \nabla f(\mathbf{x}) = \begin{bmatrix} \frac{\partial f}{\partial x_1} \\ \frac{\partial f}{\partial x_2} \\ \vdots \\ \frac{\partial f}{\partial x_n} \end{bmatrix} \]Объяснение: Градиент — это вектор, содержащий все частные производные скалярной функции. Он указывает направление наискорейшего подъема функции и широко используется в алгоритмах оптимизации МО, таких как градиентный спуск.
Пример: Для \(f(x, y) = x^2 + y^2\), \(\nabla f(x, y) = \begin{bmatrix} 2x \\ 2y \end{bmatrix}\).
Реализация (численное вычисление):
import numpy as np
def gradient(f, point, h=1e-5):
"""Вычисляет численный градиент функции f в точке point."""
point = np.asarray(point)
grad = np.zeros_like(point, dtype=float)
for i in range(len(point)):
point_plus_h = point.copy()
point_plus_h[i] += h
grad[i] = (f(*point_plus_h) - f(*point)) / h # Используем прямую разность
# Можно использовать центральную разность для лучшей точности:
# point_minus_h = point.copy()
# point_minus_h[i] -= h
# grad[i] = (f(*point_plus_h) - f(*point_minus_h)) / (2 * h)
return grad
# Пример использования:
def my_multi_func(x, y):
return x**2 + y**2
point = (1, 2)
grad_vector = gradient(my_multi_func, point)
print(f"Gradient at {point} ~ {grad_vector}") # Должно быть близко к [2, 4]
3.12. Вторая производная (Гессиан)
\[ H(f) = \begin{bmatrix} \frac{\partial^2 f}{\partial x_1^2} & \frac{\partial^2 f}{\partial x_1 \partial x_2} & \dots \\ \frac{\partial^2 f}{\partial x_2 \partial x_1} & \frac{\partial^2 f}{\partial x_2^2} & \dots \\ \vdots & \vdots & \ddots \end{bmatrix} \]Объяснение: Гессиан — это квадратная матрица вторых частных производных. Она используется в оптимизации для оценки кривизны и свойств сходимости функции.
Пример: Для \(f(x, y) = x^2 + y^2\). \(\frac{\partial^2 f}{\partial x^2} = 2\), \(\frac{\partial^2 f}{\partial y^2} = 2\), \(\frac{\partial^2 f}{\partial x \partial y} = 0\), \(\frac{\partial^2 f}{\partial y \partial x} = 0\). Гессиан \(H(f) = \begin{bmatrix} 2 & 0 \\ 0 & 2 \end{bmatrix}\).
Реализация (численное вычисление):
import numpy as np
def hessian(f, point, h=1e-5):
"""Вычисляет численный гессиан функции f в точке point."""
point = np.asarray(point)
n = len(point)
hess = np.zeros((n, n), dtype=float)
f_point = f(*point) # Значение функции в исходной точке
for i in range(n):
# Диагональные элементы (вторая производная по одной переменной)
point_plus_h = point.copy()
point_minus_h = point.copy()
point_plus_h[i] += h
point_minus_h[i] -= h
hess[i, i] = (f(*point_plus_h) - 2 * f_point + f(*point_minus_h)) / (h**2)
# Внедиагональные элементы (смешанные производные)
for j in range(i + 1, n):
point_pp = point.copy()
point_pm = point.copy()
point_mp = point.copy()
point_mm = point.copy()
point_pp[i] += h; point_pp[j] += h
point_pm[i] += h; point_pm[j] -= h
point_mp[i] -= h; point_mp[j] += h
point_mm[i] -= h; point_mm[j] -= h
hess[i, j] = (f(*point_pp) - f(*point_pm) - f(*point_mp) + f(*point_mm)) / (4 * h**2)
hess[j, i] = hess[i, j] # Гессиан симметричен для "хороших" функций
return hess
# Пример использования:
def my_multi_func(x, y):
return x**2 + y**2
point = (1, 2)
hess_matrix = hessian(my_multi_func, point)
print(f"Hessian at {point} ~ \n{hess_matrix}") # Должно быть близко к [[2, 0], [0, 2]]
3.13. Производная по направлению
\[ D_{\mathbf{v}}f(\mathbf{x}) = \nabla f(\mathbf{x}) \cdot \mathbf{v} \] (где \(\mathbf{v}\) - единичный вектор направления)Объяснение: Производная по направлению измеряет скорость изменения функции в направлении заданного вектора. Она важна в оптимизации и МО для оценки поведения функции в определенном направлении.
Пример: Для \(f(x, y) = x^2 + y^2\), \(\nabla f(x, y) = \begin{bmatrix} 2x \\ 2y \end{bmatrix}\). В направлении \(\mathbf{v} = \begin{bmatrix} 1 \\ 0 \end{bmatrix}\) (единичный вектор), \(D_{\mathbf{v}}f(x, y) = \begin{bmatrix} 2x \\ 2y \end{bmatrix} \cdot \begin{bmatrix} 1 \\ 0 \end{bmatrix} = 2x\).
Реализация:
import numpy as np
# grad_f - функция, вычисляющая градиент
def directional_derivative(grad_f, point, direction):
"""Вычисляет производную по направлению."""
grad = grad_f(point)
direction = np.asarray(direction)
# Нормализуем вектор направления, если он не единичный
norm_direction = direction / np.linalg.norm(direction)
return np.dot(grad, norm_direction)
# Пример использования (используя функцию gradient из п.3.11):
def my_multi_func(x, y):
return x**2 + y**2
point = (1, 2)
direction = [1, 0]
dir_deriv = directional_derivative(lambda p: gradient(my_multi_func, p), point, direction)
print(f"Directional derivative at {point} in direction {direction} ~ {dir_deriv}") # Близко к 2
3.14. Частные производные высших порядков
\[ \frac{\partial^k f}{\partial x_{p_1} \partial x_{p_2} \dots \partial x_{p_k}} \]Объяснение: Частные производные высших порядков расширяют понятие частных производных на порядки выше первого. Смешанные производные часто равны ( \(f_{xy} = f_{yx}\)) при условии гладкости функции (Теорема Шварца/Клеро).
Пример: Для \(f(x, y) = x^2y\). \(\frac{\partial f}{\partial y} = x^2\). \(\frac{\partial^2 f}{\partial x \partial y} = \frac{\partial}{\partial x}(x^2) = 2x\).
Реализация (численное вычисление - сложно и неточно):
Вычисление смешанных производных высших порядков численно является сложной задачей и подвержено накоплению ошибок. Часто используются библиотеки символьных вычислений (например, SymPy) или автоматического дифференцирования (например, JAX, TensorFlow, PyTorch).
# Пример для второй смешанной производной f_xy с использованием численных методов
# (аналогично вычислению внедиагональных элементов гессиана)
def mixed_partial_xy(f, point, h=1e-5):
x, y = point
f_pp = f(x + h, y + h)
f_pm = f(x + h, y - h)
f_mp = f(x - h, y + h)
f_mm = f(x - h, y - h)
return (f_pp - f_pm - f_mp + f_mm) / (4 * h**2)
# Пример использования:
def my_func(x, y):
return x**2 * y
point = (3, 4)
mixed_deriv = mixed_partial_xy(my_func, point)
print(f"d2f/dxdy at {point} ~ {mixed_deriv}") # Должно быть близко к 2*3 = 6
3.15. Полная производная
\[ \frac{df}{dt} = \sum_{i=1}^{n} \frac{\partial f}{\partial x_i} \frac{dx_i}{dt} \] (где \(f = f(x_1, \dots, x_n)\) и \(x_i = x_i(t)\))Объяснение: Полная производная учитывает изменения по всем независимым переменным, которые сами являются функциями некоторой внешней переменной \(t\). Используется в динамических системах и оптимизации.
Пример: Если \(f(x, y) = x^2 + y^2\), \(x = t\), и \(y = t^2\). Тогда \(\frac{\partial f}{\partial x} = 2x\), \(\frac{\partial f}{\partial y} = 2y\), \(\frac{dx}{dt} = 1\), \(\frac{dy}{dt} = 2t\). Полная производная \(\frac{df}{dt} = \frac{\partial f}{\partial x}\frac{dx}{dt} + \frac{\partial f}{\partial y}\frac{dy}{dt} = (2x)(1) + (2y)(2t) = 2(t) + 2(t^2)(2t) = 2t + 4t^3\).
Реализация:
# partials - список функций частных производных [df/dx1, df/dx2, ...]
# dx_dt - список функций производных [dx1/dt, dx2/dt, ...]
# point_funcs - список функций [x1(t), x2(t), ...]
# t - значение переменной t
def total_derivative(partials, dx_dt, point_funcs, t):
"""Вычисляет полную производную df/dt."""
current_point = [func(t) for func in point_funcs] # Значения x_i в момент t
total_deriv = 0
for i in range(len(partials)):
# Вычисляем частную производную в текущей точке
partial_val = partials[i](*current_point)
# Вычисляем производную dx_i/dt
dx_i_dt_val = dx_dt[i](t)
total_deriv += partial_val * dx_i_dt_val
return total_deriv
# Пример использования:
def df_dx(x, y): return 2*x
def df_dy(x, y): return 2*y
def dx_dt(t): return 1
def dy_dt(t): return 2*t
def x_func(t): return t
def y_func(t): return t**2
t_val = 3
total_deriv_val = total_derivative(
[df_dx, df_dy],
[dx_dt, dy_dt],
[x_func, y_func],
t_val
)
# 2*t_val + 4*t_val**3 = 2*3 + 4*3**3 = 6 + 4*27 = 6 + 108 = 114
print(f"Total derivative df/dt at t={t_val} is {total_deriv_val}")
3.16. Неявное дифференцирование
\[ \frac{dy}{dx} = - \frac{\partial F / \partial x}{\partial F / \partial y} \] (для уравнения \(F(x, y) = 0\))Объяснение: Неявное дифференцирование вычисляет производную зависимой переменной в уравнении, где переменная не может быть явно выражена. Используется в МО и матанализе для работы со сложными уравнениями.
Пример: Для \(F(x, y) = x^2 + y^2 - 1 = 0\) (уравнение окружности). \(\frac{\partial F}{\partial x} = 2x\), \(\frac{\partial F}{\partial y} = 2y\). Тогда \(\frac{dy}{dx} = - \frac{2x}{2y} = - \frac{x}{y}\) (при \(y \neq 0\)).
Реализация:
# F - функция F(x, y)
# partial_F_x - функция частной производной dF/dx
# partial_F_y - функция частной производной dF/dy
def implicit_differentiation(partial_F_x, partial_F_y, x, y):
"""Вычисляет dy/dx для F(x, y) = 0."""
pFy = partial_F_y(x, y)
if pFy == 0:
raise ValueError("Частная производная по y равна нулю, dy/dx не определена или бесконечна")
pFx = partial_F_x(x, y)
return -pFx / pFy
# Пример использования:
def pF_dx(x, y): return 2*x
def pF_dy(x, y): return 2*y
# Найдем производную в точке на окружности, например (sqrt(2)/2, sqrt(2)/2)
x_val = np.sqrt(2)/2
y_val = np.sqrt(2)/2
dy_dx = implicit_differentiation(pF_dx, pF_dy, x_val, y_val)
print(f"dy/dx at ({x_val}, {y_val}) is {dy_dx}") # Должно быть -x/y = -1
3.17. Разложение в ряд Тейлора
\[ f(x) \approx f(a) + f'(a)(x - a) + \frac{f''(a)}{2!}(x - a)^2 + \dots + \frac{f^{(n)}(a)}{n!}(x - a)^n + \dots \]Объяснение: Ряд Тейлора аппроксимирует функцию в окрестности точки \(a\), используя ее производные. Используется в оптимизации и численном анализе.
Пример: Для \(f(x) = e^x\) в окрестности \(a = 0\). \(f(0)=1, f'(0)=1, f''(0)=1, \dots\). Тогда \(f(x) \approx 1 + 1(x-0) + \frac{1}{2!}(x-0)^2 + \dots = 1 + x + \frac{x^2}{2} + \dots\).
Реализация (вычисление частичной суммы):
import numpy as np
import math
# derivatives - список функций производных [f(a), f'(a), f''(a), ...]
def taylor_series_approx(derivatives_at_a, a, x, terms):
"""Вычисляет приближение функции с помощью ряда Тейлора."""
result = 0
for n in range(terms):
if n >= len(derivatives_at_a): # Если не хватает посчитанных производных
print(f"Предупреждение: требуется {terms} членов, но доступно только {len(derivatives_at_a)} производных.")
break
term = derivatives_at_a[n] * (x - a)**n / math.factorial(n)
result += term
return result
# Пример использования для e^x около a=0
# Производные e^x в точке 0 равны 1
derivs_e_at_0 = [1, 1, 1, 1, 1] # f(0), f'(0), f''(0), ...
a = 0
x = 0.5
approx_e_0_5 = taylor_series_approx(derivs_e_at_0, a, x, terms=4) # Используем 4 члена (до x^3)
actual_e_0_5 = np.exp(0.5)
print(f"Taylor approx of e^0.5 (4 terms): {approx_e_0_5}")
# 1 + 0.5 + (0.5^2)/2 + (0.5^3)/6 = 1 + 0.5 + 0.125 + 0.125/6 = 1.625 + 0.02083... = 1.64583...
print(f"Actual e^0.5: {actual_e_0_5}") # ~ 1.6487
3.18. Матрица Якоби (Якобиан)
\[ J(\mathbf{f}) = \begin{bmatrix} \frac{\partial f_1}{\partial x_1} & \dots & \frac{\partial f_1}{\partial x_n} \\ \vdots & \ddots & \vdots \\ \frac{\partial f_m}{\partial x_1} & \dots & \frac{\partial f_m}{\partial x_n} \end{bmatrix} \] (для векторной функции \(\mathbf{f}: \mathbb{R}^n \to \mathbb{R}^m\))Объяснение: Матрица Якоби содержит все частные производные первого порядка векторнозначной функции. Она важна в МО для градиентной оптимизации в многомерных пространствах.
Пример: Для \(\mathbf{f}(x, y) = \begin{bmatrix} x^2 + y \\ y^2 + x \end{bmatrix}\). \(f_1 = x^2+y, f_2 = y^2+x\). \(\frac{\partial f_1}{\partial x} = 2x\), \(\frac{\partial f_1}{\partial y} = 1\). \(\frac{\partial f_2}{\partial x} = 1\), \(\frac{\partial f_2}{\partial y} = 2y\). Якобиан \(J(\mathbf{f}) = \begin{bmatrix} 2x & 1 \\ 1 & 2y \end{bmatrix}\).
Реализация (численное вычисление):
import numpy as np
# f_list - список функций [f1, f2, ..., fm]
# point - точка (x1, x2, ..., xn)
def jacobian(f_list, point, h=1e-5):
"""Вычисляет численную матрицу Якоби."""
point = np.asarray(point)
m = len(f_list)
n = len(point)
J = np.zeros((m, n), dtype=float)
for i in range(m): # по строкам (функциям fi)
f_i = f_list[i]
f_point = f_i(*point)
for j in range(n): # по столбцам (переменным xj)
point_plus_h = point.copy()
point_plus_h[j] += h
J[i, j] = (f_i(*point_plus_h) - f_point) / h # Прямая разность
# Можно использовать центральную разность
# point_minus_h = point.copy()
# point_minus_h[j] -= h
# J[i, j] = (f_i(*point_plus_h) - f_i(*point_minus_h)) / (2 * h)
return J
# Пример использования:
def f1(x, y): return x**2 + y
def f2(x, y): return y**2 + x
point = (1, 2)
jac_matrix = jacobian([f1, f2], point)
print(f"Jacobian at {point} ~ \n{jac_matrix}") # Должно быть близко к [[2*1, 1], [1, 2*2]] = [[2, 1], [1, 4]]
3.19. Длина дуги кривой
\[ L = \int_{a}^{b} \sqrt{1 + \left( \frac{dy}{dx} \right)^2} dx \] (для кривой \(y = f(x)\))Объяснение: Длина дуги измеряет расстояние вдоль кривой между двумя точками. Используется в геометрии и физике для анализа пути.
Пример: Для \(y = x^2\) на отрезке \([0, 1]\). \(\frac{dy}{dx} = 2x\). \(L = \int_{0}^{1} \sqrt{1 + (2x)^2} dx = \int_{0}^{1} \sqrt{1 + 4x^2} dx\).
Реализация (численное интегрирование):
from scipy.integrate import quad
import numpy as np
# f_prime - функция, вычисляющая производную dy/dx
def arc_length(f_prime, a, b):
"""Вычисляет длину дуги кривой y=f(x) от a до b."""
integrand = lambda x: np.sqrt(1 + f_prime(x)**2)
length, error_estimate = quad(integrand, a, b)
return length
# Пример использования:
def y_prime(x):
return 2*x
length_val = arc_length(y_prime, 0, 1)
print(f"Arc length of y=x^2 from 0 to 1 is ~ {length_val}") # ~ 1.4789
3.20. Кривизна функции
\[ \kappa(x) = \frac{|y''(x)|}{(1 + [y'(x)]^2)^{3/2}} \] (для кривой \(y = f(x)\))Объяснение: Кривизна количественно определяет, насколько резко изгибается кривая в данной точке. Используется в геометрии и анализе траекторий в робототехнике и МО.
Пример: Для \(y = x^2\). \(y'(x) = 2x\), \(y''(x) = 2\). Тогда \(\kappa(x) = \frac{|2|}{(1 + (2x)^2)^{3/2}} = \frac{2}{(1 + 4x^2)^{3/2}}\).
Реализация:
import numpy as np
# f_prime - функция производной y'
# f_double_prime - функция второй производной y''
def curvature(f_prime, f_double_prime, x):
"""Вычисляет кривизну функции y=f(x) в точке x."""
y_prime_val = f_prime(x)
y_double_prime_val = f_double_prime(x)
numerator = abs(y_double_prime_val)
denominator = (1 + y_prime_val**2)**1.5
if denominator == 0:
# Это может произойти, если y'(x) - комплексное число или в особых точках,
# но для реальных функций знаменатель >= 1.
# Возвращаем бесконечность или NaN, или обрабатываем иначе.
return np.inf
return numerator / denominator
# Пример использования:
def y_prime(x): return 2*x
def y_double_prime(x): return 2 # Можно сделать функцией: lambda x: 2
curvature_at_0 = curvature(y_prime, y_double_prime, 0)
print(f"Curvature of y=x^2 at x=0 is {curvature_at_0}") # 2 / (1+0)^1.5 = 2
curvature_at_1 = curvature(y_prime, y_double_prime, 1)
print(f"Curvature of y=x^2 at x=1 is {curvature_at_1}") # 2 / (1+4)^1.5 = 2 / 5^1.5 ~ 0.1789
3.21. Интегрирование по частям
\[ \int u \, dv = uv - \int v \, du \] (где \(dv = v'(x)dx\), \(du = u'(x)dx\))Объяснение: Интегрирование по частям — это техника, производная от правила произведения для дифференцирования. Она используется для упрощения интегралов, включающих произведения функций.
Пример: Для \(\int x e^x dx\). Пусть \(u = x\) и \(dv = e^x dx\). Тогда \(du = dx\) и \(v = \int e^x dx = e^x\). Применяя формулу: \(\int x e^x dx = x e^x - \int e^x dx = x e^x - e^x + C\).
Реализация (символьное интегрирование):
import sympy
x_sym = sympy.symbols('x')
u = x_sym
f = sympy.exp(x_sym) # Полная функция под интегралом f = x * exp(x)
# SymPy может интегрировать по частям автоматически
integral_expr = sympy.integrate(u * f, x_sym) # Sympy попытается найти интеграл
print(integral_expr) # Должен вывести (x - 1)*exp(x)
# Или явно указать u и dv для функции integrate_by_parts (если она доступна)
# В стандартном sympy может не быть отдельной integrate_by_parts,
# но можно эмулировать шаги:
# dv = sympy.exp(x) * dx -> v = sympy.integrate(sympy.exp(x), x) = exp(x)
# du = sympy.diff(x, x) * dx = 1 * dx
# integral = u * v - sympy.integrate(v * sympy.diff(u, x), x)
# print(integral)
Примечание: Как и в случае с интегралом степенной функции, этот пример показывает правило символьно.
3.22. Объем тела вращения (Метод дисков)
\[ V = \pi \int_{a}^{b} [f(x)]^2 dx \] (при вращении кривой \(y=f(x)\) вокруг оси Ox)Объяснение: Метод дисков вычисляет объем тела, полученного вращением плоской фигуры вокруг оси, путем разбиения его на тонкие диски, перпендикулярные оси вращения. Распространен в геометрии и физике.
Пример: Для \(y = x^2\), вращаемой вокруг оси x на отрезке \([0, 1]\). \(V = \pi \int_{0}^{1} (x^2)^2 dx = \pi \int_{0}^{1} x^4 dx = \pi \left[ \frac{x^5}{5} \right]_{0}^{1} = \pi (\frac{1^5}{5} - \frac{0^5}{5}) = \frac{\pi}{5}\).
Реализация (численное интегрирование):
from scipy.integrate import quad
import numpy as np
def volume_of_revolution(f, a, b):
"""Вычисляет объем тела вращения (вокруг оси x) методом дисков."""
integrand = lambda x: np.pi * f(x)**2
volume, error_estimate = quad(integrand, a, b)
return volume
# Пример использования:
def my_func(x):
return x**2
volume_val = volume_of_revolution(my_func, 0, 1)
print(f"Volume of revolution for y=x^2 on [0,1] is ~ {volume_val}") # Близко к pi/5 ~ 0.628
print(f"Pi/5 = {np.pi/5}")
3.23. Поверхностный интеграл (1-го рода)
\[ \iint_S f(x, y, z) dS = \iint_R f(x, y, g(x, y)) \sqrt{1 + \left(\frac{\partial g}{\partial x}\right)^2 + \left(\frac{\partial g}{\partial y}\right)^2} dA \] (где поверхность S задана уравнением \(z = g(x, y)\), R - проекция S на плоскость xy)Объяснение: Поверхностный интеграл расширяет идею криволинейного интеграла на поверхность, суммируя скалярное поле или поток векторного поля по поверхности.
Пример: Вычислить поверхностный интеграл \(f(x, y, z) = z\) по поверхности \(z = x^2 + y^2\) для \(x^2 + y^2 \le 1\). Здесь \(g(x,y) = x^2+y^2\). \(\frac{\partial g}{\partial x} = 2x\), \(\frac{\partial g}{\partial y} = 2y\). Интеграл равен \(\iint_R (x^2+y^2) \sqrt{1+(2x)^2+(2y)^2} dA\), где R - круг \(x^2+y^2 \le 1\). Проще вычислить в полярных координатах.
Реализация (численное интегрирование):
from scipy.integrate import dblquad
import numpy as np
# f - скалярная функция f(x, y, z)
# g - функция поверхности z=g(x, y)
# g_dx - функция частной производной dg/dx
# g_dy - функция частной производной dg/dy
# bounds_x - функция или кортеж пределов для x (могут зависеть от y)
# bounds_y - кортеж пределов для y
def surface_integral_scalar(f, g, g_dx, g_dy, bounds_y, bounds_x):
"""Вычисляет поверхностный интеграл 1-го рода."""
def integrand(y, x): # Порядок аргументов важен для dblquad (сначала внутренняя переменная)
z = g(x, y)
gx = g_dx(x, y)
gy = g_dy(x, y)
surface_element = np.sqrt(1 + gx**2 + gy**2)
return f(x, y, z) * surface_element
# dblquad(func, a, b, gfun, hfun) вычисляет int_a^b dx int_{gfun(x)}^{hfun(x)} dy func(y, x)
# Нам нужно int dy int dx integrand(y, x)
# Если пределы const: dblquad(integrand, y_min, y_max, lambda y: x_min, lambda y: x_max) - НЕПРАВИЛЬНО
# Если пределы const: dblquad(integrand, x_min, x_max, lambda x: y_min, lambda x: y_max) - ПРАВИЛЬНО
# bounds_y = (y_min, y_max), bounds_x = (x_min, x_max) или функции
if callable(bounds_x[0]): # Если пределы x зависят от y
integral_val, error = dblquad(integrand, bounds_y[0], bounds_y[1], bounds_x[0], bounds_x[1])
else: # Если пределы константы
integral_val, error = dblquad(integrand, bounds_y[0], bounds_y[1], lambda y: bounds_x[0], lambda y: bounds_x[1])
# ИЛИ, поменяв порядок интегрирования, если удобнее:
# def integrand_yx(x, y): ...
# integral_val, error = dblquad(integrand_yx, bounds_x[0], bounds_x[1], lambda x: bounds_y[0], lambda x: bounds_y[1])
return integral_val
# Пример использования для z по z=x^2+y^2 над кругом x^2+y^2 <= 1
def f(x, y, z): return z
def g(x, y): return x**2 + y**2
def g_dx(x, y): return 2*x
def g_dy(x, y): return 2*y
# Пределы для круга x^2+y^2 <= 1: y от -1 до 1, x от -sqrt(1-y^2) до sqrt(1-y^2)
bounds_y = (-1, 1)
bounds_x = (lambda y: -np.sqrt(1 - y**2), lambda y: np.sqrt(1 - y**2))
# Вычисление может быть сложным из-за пределов и функции
try:
surf_integral_val = surface_integral_scalar(f, g, g_dx, g_dy, bounds_y, bounds_x)
print(f"Surface integral value ~ {surf_integral_val}")
except Exception as e:
print(f"Could not compute integral numerically: {e}")
# Переход в полярные координаты часто упрощает вычисление:
# x=r*cos(theta), y=r*sin(theta), z=r^2
# dA = r dr dtheta
# sqrt(1+4x^2+4y^2) = sqrt(1+4r^2)
# Integral = int_0^{2pi} dtheta int_0^1 (r^2) * sqrt(1+4r^2) * r dr
integrand_polar = lambda r: r**3 * np.sqrt(1 + 4*r**2)
integral_r, _ = quad(integrand_polar, 0, 1)
surf_integral_val_polar = 2 * np.pi * integral_r
print(f"Surface integral value (polar) ~ {surf_integral_val_polar}") # ~ 3.829
3.24. Дивергенция векторного поля
\[ \text{div } \mathbf{F} = \nabla \cdot \mathbf{F} = \frac{\partial F_1}{\partial x} + \frac{\partial F_2}{\partial y} + \frac{\partial F_3}{\partial z} \] (где \(\mathbf{F} = (F_1, F_2, F_3)\))Объяснение: Дивергенция измеряет величину источника или стока векторного поля в данной точке. Используется в гидродинамике и электромагнетизме.
Пример: Для \(\mathbf{F} = \begin{bmatrix} x \\ y \\ z \end{bmatrix}\). \(\text{div } \mathbf{F} = \frac{\partial x}{\partial x} + \frac{\partial y}{\partial y} + \frac{\partial z}{\partial z} = 1 + 1 + 1 = 3\).
Реализация (символьное вычисление):
import sympy
x, y, z = sympy.symbols('x y z')
F = sympy.Matrix([x, y, z]) # Векторное поле как матрица sympy
# Вычисляем дивергенцию
divergence = sympy.diff(F[0], x) + sympy.diff(F[1], y) + sympy.diff(F[2], z)
print(f"Divergence of F = [x, y, z] is: {divergence}") # Выведет 3
# Для другого поля
F2 = sympy.Matrix([x*y, y*z, z*x])
div2 = sympy.diff(F2[0], x) + sympy.diff(F2[1], y) + sympy.diff(F2[2], z)
print(f"Divergence of F2 = [xy, yz, zx] is: {div2}") # Выведет y + z + x
Примечание: Численное вычисление дивергенции потребовало бы численного вычисления частных производных.
3.25. Ротор (вихрь) векторного поля
\[ \text{curl } \mathbf{F} = \nabla \times \mathbf{F} = \begin{vmatrix} \mathbf{i} & \mathbf{j} & \mathbf{k} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ F_1 & F_2 & F_3 \end{vmatrix} \]Объяснение: Ротор измеряет вращение или циркуляцию векторного поля в точке. Критически важен в гидромеханике и электромагнетизме.
Пример: Для \(\mathbf{F} = \begin{bmatrix} 0 \\ 0 \\ xy \end{bmatrix}\). \(\text{curl } \mathbf{F} = \mathbf{i} (\frac{\partial(xy)}{\partial y} - \frac{\partial 0}{\partial z}) - \mathbf{j} (\frac{\partial(xy)}{\partial x} - \frac{\partial 0}{\partial z}) + \mathbf{k} (\frac{\partial 0}{\partial x} - \frac{\partial 0}{\partial y})\) \( = \mathbf{i}(x - 0) - \mathbf{j}(y - 0) + \mathbf{k}(0 - 0) = x\mathbf{i} - y\mathbf{j} = \begin{bmatrix} x \\ -y \\ 0 \end{bmatrix}\). (Примечание: в исходном OCR примере ответ был [-y, x, 0], что соответствует полю F = [-y, x, 0] или F=[0,0,z(x^2+y^2)/2]. Для F=[0,0,xy] ответ [x, -y, 0]. Уточним по OCR реализации. OCR реализация использует `F.jacobian().transpose() - F.jacobian()`, что не является стандартной формулой ротора.)
Реализация (символьное вычисление):
import sympy
x, y, z = sympy.symbols('x y z')
F1, F2, F3 = 0, 0, x*y # Компоненты поля F = [0, 0, xy]
# Вычисляем компоненты ротора
curl_x = sympy.diff(F3, y) - sympy.diff(F2, z)
curl_y = sympy.diff(F1, z) - sympy.diff(F3, x)
curl_z = sympy.diff(F2, x) - sympy.diff(F1, y)
curl_F = sympy.Matrix([curl_x, curl_y, curl_z])
print(f"Curl of F = [0, 0, xy] is: {curl_F}") # Выведет Matrix([[x], [-y], [0]])
# Используя встроенную функцию sympy (если доступна в вашей версии)
# from sympy.vector import curl, CoordSys3D
# N = CoordSys3D('N')
# F_vec = F1*N.i + F2*N.j + F3*N.k
# curl_sympy = curl(F_vec)
# print(f"Curl using sympy.vector: {curl_sympy}")
Примечание: Реализация из OCR использует разность якобиана и его транспонирования, что связано с вихрем, но не является прямым вычислением ротора по стандартной формуле определителя. Выше приведена реализация по стандартной формуле.