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


Введение

Математический анализ знает множество прекрасных функций с необычными свойствами. Среди них интегральные синус 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, а приx> e берем 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-функции Ламберта могут быть достаточно точно вычислены за короткое время, что делает возможным решение некоторых типов уравнений.




Комментарии (11):

  1. KvanTTT
    /#24341742 / +7

    В детстве в школе "изобрел" тетрацию или гипероператор-4 (только потом узнал как это называется). И только когда вырос, узнал, что некоторые уравнения с тетрацией решаются с использованием W-функции ламберта. Самое простое:


    x^x = x, x - ?
    x = ln(k) / W(ln(k))

    Так что эта функция вызывает у меня приятные ностальгические воспоминания :)


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

    Разве не более эффективней раскладывать в ряд Тейлора и не использовать ни логарифмы, ни рекурсию? https://en.wikipedia.org/wiki/Lambert_W_function#Asymptotic_expansions

    • Sdima1357
      /#24342010

      x^x=x для х не равного 0 имеем х^(х-1) = 1 , имеет два вполне очевидных корня +1 и -1 или я что-то путаю?

      • KvanTTT
        /#24342036 / +1

        Извините, опечатался. Имел в виду x ^ x = k, где k — константа. В решении написано правильно с k.

    • e1vanov
      /#24342290

      Доброй ночи! Я использую рекурсию и логарифмы при подсчете значений так называемой -1-ветви функции Ламберта, а разложение в ряд Тейлора дает значение 0-ветви в окрестности нуля

      • KvanTTT
        /#24342950

        Кстати, можно переписать функцию так, чтобы не использовать рекурсию, так эффективней, как-то так:


        def LambertW_1(x):
            result = log(-x)
            for t in range(99):
                result = log(x / result)
            return result

    • KvanTTT
      /#24342930

      Если использовать тетрацию, то запись уравнения выглядела бы так: x^^2 = k.

    • alxndrlsn
      /#24344408 / +1

      любопытно... а меня с некоторых пор "тревожила" идея "мультипликата", в основе которой лежит предположение о том, что если для интегрируемой функции g(x) имеется ее первообразная G(x), то для функции f(x) представленной в виде экспоненты в степени g(x), тоже должна быть некая первообразная F(x) (и может даже только одна?), получаемая из f(x) некоторым преобразованием ("мультиплицированием"), которое соответствует интегрированию степени экспоненты (прошу прощения за невнятное пояснение, лучше взгляните на рисунок - осторожно! может быть заразно!):

      M и ? - это и есть те недоступные моему пониманию функции
      M и ? - это и есть те недоступные моему пониманию функции

      Может быть, кто-то тоже переболел такими "фантазиями" и разрешился от них?

  2. belch84
    /#24342648 / +2

    Можно подойти к решению уравнения
    image
    немного иначе.
    Функция Ламберта является решением обыкновенного дифференциального уравнения (ОДУ)
    image
    Если доступно средство для численного решения дифференциальных уравнений, можно найти его решение, и вычислить корень по приведенной в публикации формуле
    image
    Конечно, численное решение дифференциального уравнения устроено намного сложнее, и при этом трудно достигнуть очень высокой точности, но при наличии готовой программы для решения ОДУ можно решить исходное уравнение совсем без программирования. Для уравнения
    image
    у меня получился корень 1.03215

    Сравнение точного графика ф-ции Ламберта и решения дифференциального уравнения
    image
    Зеленый — график функции Ламберта, красный — решение дифференциального уравнения

  3. Un_ka
    /#24342860

    Графики были построены в GeoGebra?

    • e1vanov
      /#24343548

      В том числе. Еще использовались desmos и встроенное в macOS приложение "графический калькулятор"