W-функция Ламберта и ее приложения

Введение

Математический анализ знает множество прекрасных функций с необычными свойствами. Среди них интегральные синус si(x)и логарифм li(x), также нельзя не отметить гамма-функцию \Gamma(x) или очень известную дзета-функцию Римана \zeta(s). Но сегодня я предлагаю читателю посмотреть на функцию W-Ламберта W(x).

Что такое W-функция Ламберта?

Для того, чтобы понять, что такое W-функция Ламберта, достаточно посмотреть на следующее равенство, которое по аналогии с основным тригометрическим тождеством предлагаю именовать «основное Ламбертово тождество»:

W(x)e^{W(x)} = x

Другими словами, функция Ламберта — обратная для f(x) =xe^x. Однако после первых же исследований станет понятно, что f не инъективна, а именно одно и то же значение y достигается при двух разных аргументах, если y \in (-e^{-1}; 0). Поэтому данное выше определение требует пояснений.

График y = xe^x
График y = xe^x

Исследовав производную функции f' = e^x(x+1), понимаем, что функция возрастает на [-1; +\infty)и убывает на (-\infty; -1]. Таким образом, давайте построим обратную функцию к данной на соответствующих промежутках монотонности.

Графики W-функции Ламберта (черный) и y = xe^x (серый)
Графики W-функции Ламберта (черный) и y = xe^x (серый)

Ветвь, для которой y\in [-1; +\infty), называется W_{0}(x), другая — W_{-1}(x).

Постановка задачи

Задача. Научиться находить действительные корни уравнения следующего вида:

x^pe^x=q \\p \in \mathbb{Q} \ \backslash \ \mathbb{Z}, q\in\mathbb{R_{+}}

Как вы, наверное, уже догадались, для решения будем использовать W-функцию Ламберта. Итак, сначала возведем и левую, и правую часть в степень p^{-1}(это преобразование не является равносильным при четных целых p , а при нечетных нужно будет расширять множество значений q до всех действительных чисел, поэтому решаем задачу для указанных выше ограничений).

xe^{\frac{x}p} = q^{\frac{1}{p}}

Теперь для того, чтобы воспользоваться основным Ламбертовым тождеством, нам нужно получить выражение с x такое, как и в показателе степени экспоненты. Для этого поделим и левую, и правую часть на p.

\frac{x}{p}e^{\frac{x}{p}}=\frac{q^{\frac{1}{p}}}{p}

И теперь мы можем воспользоваться основным Ламбертовым тождеством:

\frac{x}{p}=W(\frac{q^\frac{1}{p}}{p})

Откуда и получаем уже итоговую формулу для x.

x = p*W(\frac{q^\frac{1}{p}}{p})

Вычисление W-функции Ламберта

Заметим, что при x\in (-e^{-1}; 0) функция Ламберта дает два действительных значения: по одному на каждой из ветвей W_{0}(x) и W_{-1}(x) соответственно. В этом случае у изначального уравнения будет 2 корня.

Вычисление W0. Будем использовать метод бинарного поиска по ответу. Мы можем так поступить, поскольку W_{0}(x) возрастает на [-e^{-1}; +\infty).

Левая граница бинарного поиска понятна и равна-1. Теперь возникает вопрос, как выбрать правую границу. Первая идея, которая приходит на ум: положить ее равной x, ибо верно неравенство x \geq W(x), причем равенство достигается только в нуле.

x \geq W(x) \iff x \geq \frac{x}{e^{W(x)}} \iff (e^{W(x)}-1)x \geq0

Однако для достаточно больших x это может быть не лучшим вариантом. Поэтому давайте посмотрим на другой: выберем правую границу равной ln(x).

ln(x) \geq W(x) \iff x \geq e^{W(x)} \iff x\geq \frac{x}{W(x)} \iff x(\frac{W(x)-1}{W(x)})\geq0 \\ x\geq e

Итого: при-e^{-1}\leq x\leq e правой границей выбираем x, а приln(x).

Графики y = x (синий), y = W(x) (зеленый) и y = ln(x) (красный)
Графики y = x (синий), y = W(x) (зеленый) и y = ln(x) (красный)

Ассимптотика: O(\frac{ln(x)}{prec}), prec — изначально выбранная точность (например, 10-12)

Вычисление W-1. Здесь будем использовать следующее бесконечное выражение для W_{-1}(x):

W_{-1}(x)=ln\frac{-x}{-ln\frac{-x}{\dots}}

Чем глубже мы спускаемся, тем выше точность вычислений.

Реализация на Python

from math import *  def LambertW0(x):     left = -1     right = x if x <= e else log(x)     prec = 10**-12 # точность     # бинарный поиск     while right - left > prec:         mid = (right + left) / 2         if mid * exp(mid) > x:             right = mid         else:             left = mid     return right  def LambertW_1(x, t): # t - показатель точности     if t == 100:       return log(-x)     else:       return log((-x)/(-LambertW_1(x, t + 1)))  def sol(p, q):     s = q**(1/p) / p     if s < -exp(-1):         return "No real solutions"     ans = "Solutions: " + str(p * LambertW0(s)) + " "     if -exp(-1) < s and s < 0:       ans += str(p * LambertW_1(s, 0))     return ans       p = float(input()) q = float(input()) print(sol(p, q))

Проверка

Уравнение 1. x^{-\frac{3}{2}}e^x = 2

Графики y = x^(-1.5)*e^x (фиолетовый) и y = 2 (синий)
Графики y = x^(-1.5)*e^x (фиолетовый) и y = 2 (синий)
Вывод программы
Вывод программы

Уравнение 2. x^{\frac{21}{10}}e^x=3

Графики y = x^2.1*e^x (фиолетовый) и y = 3 (синий)
Графики y = x^2.1*e^x (фиолетовый) и y = 3 (синий)
Вывод программы
Вывод программы

Уравнение 3. x^{-\frac{3}{4}}e^x=4

Графики y = x^(-0.75)*e^x (фиолетовый) и y = 4 (синий)
Графики y = x^(-0.75)*e^x (фиолетовый) и y = 4 (синий)
Вывод программы
Вывод программы

Заключение / выводы

Значения на 0 и -1 ветви W-функции Ламберта могут быть достаточно точно вычислены за короткое время, что делает возможным решение некоторых типов уравнений.


ссылка на оригинал статьи https://habr.com/ru/post/665634/

Добавить комментарий

Ваш адрес email не будет опубликован.