Обратная функция

[B@R5uk]

Заинтересовала меня тут одна функция. Она является решением вот этого уравнения:

$$\exp\left(y\right)+y=x$$

Её график (чёрным):


Она практически линейна при отрицательных значениях аргумента и очень похожа на логарифм при положительных. Поэтому я обзываю её линейный логарифм и обозначаю $\operatorname{lilg}\left(x\right)$. Эта функция позволяет решать уравнение

$$\exp\left(x\right)+kx=b$$

которое возникает при использовании формулы Шокли для честного обсчёта схем с диодом. Решением будет:

$$x=\ln k+\operatorname{lilg}\left(\frac{b}{k}-\ln k\right)$$

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

$$\operatorname{W}\left( x \right)=\ln x-\operatorname{lilg}\left( \ln x \right)$$

$$\operatorname{lilg}\left( y \right)=y-\operatorname{W}\left( \exp \left( y \right) \right)$$

Однако, для решения уравнения выше через W-функцию Ламберта нужно сначала брать экспоненту от аргумента, а потом вычислять логарифмоподобную W-функцию от этой экспоненты. Если аргумент будет слишком большой, то наступит переполнение числа с плавающей точкой. Ну и линейный логарифм ведёт себя хорошо на всей действительной оси. Плюс в MATLABе нет W-функции Ламберта, поэтому чем реализовывать проблематичную функцию, лучше сразу реализовать ту, что нужно.

Для расчёта значений функции я не смог придумать ничего лучше следующего. Взять начальное приближение

$${{y}_{0}}=\left\{ \begin{array}{{l}{l}} x-\exp \left( x \right), & x\in \left( -\infty ,-2 \right) \\ P\left( x \right), & x\in \left[ -2,6 \right] \\ \left( 1-\frac{1}{x} \right)\ln \left( x \right), & x\in \left( 6,+\infty \right) \\ \end{array} \right.$$

где $P\left( x \right)=0.0037{{x}^{3}}-0.0698{{x}^{2}}+0.6300x-0.5682$ и затем применить три итерации метода Ньютона:

$${{y}_{n+1}}=\frac{\left( {{y}_{n}}-1 \right)\exp \left( {{y}_{n}} \right)+x}{\exp \left( {{y}_{n}} \right)+1}$$

После этого значения приближений в пределах ошибок округления находятся вблизи реальных значений функции. Однако, "шум" ошибок округления при больших положительных значениях аргумента значительно больше того, чего бы мне хотелось.

Производная и интеграл обратной функции:

$$\operatorname{lil}{g}'\left( x \right)=\frac{1}{\exp \left( \operatorname{lilg}\left( x \right) \right)+1}$$

$$\int{\operatorname{lilg}\left( x \right)dx}=-\frac{1}{2}{{\operatorname{lilg}}^{2}}\left( x \right)+\left( x+1 \right)\operatorname{lilg}\left( x \right)-x+C$$

При изучении значения функции обычно её раскладывают в ряд Тейлора в какой-нибудь хорошей точке. Я тоже так попытался сделать:

$$\operatorname{lilg}\left( 1+x \right)=\frac{x}{2\cdot 1!}-\frac{{{x}^{2}}}{{{2}^{3}}\cdot 2!}+\frac{{{x}^{3}}}{{{2}^{5}}\cdot 3!}+\frac{{{x}^{4}}}{{{2}^{7}}\cdot 4!}-\frac{13{{x}^{5}}}{{{2}^{9}}\cdot 5!}+\frac{47{{x}^{6}}}{{{2}^{11}}\cdot 6!}+\ldots $$

Однако дальше этот ряд продлевать весьма проблематично. Формула обращения ряда настоящая жесть, которая вряд ли применима для практического счёта, а расчёт производных в лоб даёт нарастающие по сложности полиномы:

$${{\operatorname{lilg}}^{\left( n \right)}}\left( x \right)=\frac{{{P}_{n}}\left( {{e}^{y}} \right)}{{{\left( {{e}^{y}}+1 \right)}^{2n-1}}}$$

$${{P}_{1}}\left( z \right)=1$$

$${{P}_{n}}\left( z \right)=-\left( 2n-3 \right)z{{P}_{n-1}}\left( z \right)+\left( {{z}^{2}}+z \right){{{P}'}_{n-1}}\left( z \right)$$

Усмотреть закономерность в коэффициентах этих полиномов пока не получилось. Плюс, если даже получится, надо будет эти коэффициенты просуммировать для $z=1$, чтобы получить ответ в общем виде.

Однако уже с этими членами ряда можно найти значение функции в нуле с точностью 5 значащих цифр. Это постоянная Омега с обратным знаком (что бы эта постоянная не обозначала):

$$\operatorname{W}\left( 1 \right)=-\operatorname{lilg}\left( 0 \right)=\Omega =0.56714329040978\ldots $$

Казалось бы сходимость ряда очень быстрая, но не такая быстрая как у экспоненты (хотя гораздо лучше, чем у логарифма). Об ограниченности круга сходимости можно догадаться, если задуматься, однолистна ли функция в комплексной плоскости. Если по-внимательней рассмотреть уравнение $\exp(a+i b)+a+i b=x+i y$, можно понять, что это не так. То есть функция должна иметь где-то на комплексной плоскости особые точки — точки ветвления. Я попытался просчитать значения функции в комплексной плоскости численно с помощью метода Ньютона, двигаясь по чуть-чуть от действительной оси в сторону положительного направления мнимой оси. Начальными приближениями в методе Ньютона для следующего значения мнимой части аргументов становились значения функции для предыдущего значения мнимой части аргумента. Получились такие картинки:

(Действительная часть)

(Мнимая часть)

Очень похоже на то, что точками ветвления являются точки $-1\pm i\pi $.

В связи с проделанным возникает куча вопросов. Например, не изучена ли эта функция кем-нибудь давным давно и не дано ли ней какое-нибудь название? Я понимаю, что уже есть W-функция Ламберта-Эйлера, и не стоит плодить сущности без необходимости, но всё же. Может быть и разложение в ряд Тейлора в общем виде есть? Как вообще рассчитывают значения специальной функции в произвольной точке комплексной плоскости? Как строго доказать положение точек ветвления?

[Singular]

Посмотрите французскую версию Wiki, там рассмотрен Ваш случай.

[StaticZero]

Если убрать строгость, то рецепт такой. Для функции $f(x) = x - y(x)$ имеет место уравнение

$$ f'(1 + f) = f, \quad f(0) = \Omega $$

и можно получить представление

$$ f(x) \exp f(x) = \Omega e^{x+\Omega}. $$

Теперь вроде понятно, что $f \to 0$ с экспоненциальной скоростью при $x \to -\infty$. Разложим в ряд левую часть в окрестности точки $f = 0$ как $f$, получим нулевое приближение

$$ f_0(x) = \Omega e^{x+\Omega}. $$

Представим теперь $f_1= f_0 + \varepsilon$. Справа у нас нулевое приближение. Получаем

$$ f_0 = e^{f_0 + \varepsilon} (\varepsilon + f_0) \approx e^{f_0} f_1, $$

откуда $f_1 = f_0 \exp \left(- f_0 \right)$. Следующее приближение $f_2 = f_1 + \varepsilon$ строим аналогично:

$$ f_0 = e^{f_1 + \varepsilon} (\varepsilon + f_1) \approx e^{f_1} f_2, $$

общая процедура

$$ f_{n+1} = \Omega e^{x+\Omega} \exp \left( - f_n \right), $$

итерации повторять до полного просветления. Для $x \leqslant 0$ итерации точно сойдутся к $f(x)$, для $x \geqslant 0$ не знаю, может быть, вы увидите.

[B@R5uk]

А верно, что точки ветвления обратной функции связаны с точками, в которых производная прямой функции обращается в ноль? Например, синус и косинус равны $\pm 1$ в точках, где их производные равны нулю, и эти точки $\pm 1$ как раз являются точками ветвления обратных к синусу и косинусу функций. С экспонентой и логарифмом, правда, ситуация совсем другая...

-- 22.01.2019, 22:08 --

StaticZero, что-то совсем не понял, что вы решаете.

-- 22.01.2019, 22:29 --

Просчитал тут чуть-чуть дальше в глубь комплексной плоскости, судя по всему, точки $-1\pm i\pi$ так же являются точками ветвления. Только либо я что-то не правильно считаю, либо ветвление происходит не на всех листах многолистной функции. Однако это подтверждает мою гипотезу. Прямая функция: $y(x)=\exp(x)+x$, её производная: $y'(x)=\exp(x)+1$. Производная обращается в нуль в точках $x_k=\pm(2k-1)\pi$, $k=1,2,...$. В этих точках прямая функция принимает значения $y_k=-1\pm(2k-1)\pi$, то есть как раз те, что у меня получаются из численного счёта.

[StaticZero]

B@R5uk в сообщении #1370887 писал(а):

вот этого уравнения:

$$\exp\left(y\right)+y=x$$

Я вот это уравнение решаю, строю численную процедуру (вы, кажется, жаловались, что ваши вычисления неточны). А вы?

-- 22.01.2019 в 22:36 --

StaticZero в сообщении #1370954 писал(а):

для $x \geqslant 0$ не знаю, может быть, вы увидите.

Я пока не вижу, что могло бы мешать. Сходиться, может, и будут, но очень медленно.

Есть ещё один минус, значение $\Omega$ надо откуда-то брать, то есть процедуру нельзя применить для вычисления собственного параметра. Но если его задать каким надо, то для $x \leqslant 0$ сходиться действительно будет сверхбыстро.

[StaticZero]

Процесс

$$ y_{n+1} = \ln \frac{x}{1 + y_n \exp(-y_n)} $$

даёт решение для $x \geqslant 1$. Сходимость по моим наблюдениям тоже очень быстрая.

Для интервала, содержащего $[0, 1]$, подойдёт метод $y_{n+1} = x - \exp(y_n)$, $y_0 = \frac{x - 1}{2}$.

[B@R5uk]

StaticZero в сообщении #1370975 писал(а):

вы, кажется, жаловались, что ваши вычисления неточны

Не совсем так. Метод Ньютона сходится к решению и делает это ошеломляюще быстро: всего три итерации. Тут больше вопрос как численно проверить точность решения в условиях ограниченной точности компьютерной арифметики. Поскольку точного решения я не знаю, то я подставлял свой результат $\hat{y}$ в выражение:

$$err=\exp(y)+y-x$$

И сравнивал результат со значением

$$\operatorname{eps}(y)$$

где $\operatorname{eps}$ — функция, выдающая абсолютное значение последнего значащего бита мантиссы. Обычно это $2^{-52}$ помножить на модуль аргумента. Результат этой проверки для больших положительных $x$ оказывался сомнительным:



Но сейчас, я, кажется, начинаю понимать почему. Всё дело в экспоненте.

Пусть $\hat{y}=y(1+\varepsilon)$, где $y$ — точное решение, а $\varepsilon$ — относительная погрешность приближения $\hat{y}$, тогда:

$$\exp(\hat{y})=\exp(y)\exp(1+\varepsilon)\approx\exp(y)(1+\varepsilon)$$

Вроде бы относительная погрешность при возведении в степень не меняется, но при моей проверке эта экспонента завышает ошибку:

$$err=\exp(\hat{y})+\hat{y}-x=\exp(y)+y-x+\exp(y)\varepsilon+y\varepsilon=(\exp(y)+y)\varepsilon$$

Поскольку при больших положительных $x$ выполняется $\exp(y)\gg y$, то я фактически сравниваю $\exp(y)\varepsilon$ (синие точки) с $y\varepsilon$ (красная линия). Не удивительно, что получается завышенных шум ошибок округления.

Если для больших положительных $x$ использовать другую проверку, то всё получается значительно более прилично:



Для такой проверки можно записать следующее:

$$err=\ln \left( x-\widehat{y} \right)-\widehat{y}=\ln \left( x-y \right)-y+\ln \left( 1-\frac{\varepsilon }{x-y} \right)-\varepsilon \approx -\frac{\varepsilon }{x-y}-\varepsilon \approx -\varepsilon $$

То есть лучшего и желать не надо.

Меня сейчас интересует такой вопрос: как вычислять значение функции в произвольной точке комплексной плоскости? Метод Ньютона хорош, и функция для него подходящая, но он требует угадать начальное приближение с точностью порядка $0.01$, чтобы он сошёлся к решению (с точностью $2^{-52}$) за три итерации. Можно ли сделать это для всей комплексной плоскости?

[StaticZero]

B@R5uk в сообщении #1371713 писал(а):

То есть лучшего и желать не надо.

Хорошо. Если у вас относительная погрешность $\varepsilon$, то абсолютная будет $x\varepsilon$, об этом только нужно помнить.

Касательно $\mathbb C$. У меня нет матлаба под рукой проверить. А что если просто в тупую взять процедуры, которые я привёл? Вдруг сработают... (я их реализовал на Си и удовлетворился тем, что сходятся при начальном приближении, зависящем только от $x$, но в $\mathbb R$). А если сработают не везде, то тогда будет повод задуматься сильнее, что сделать нужно.

[B@R5uk]

StaticZero в сообщении #1371720 писал(а):

об этом только нужно помнить.

Да, это я нехорошо сделал. В случае логарифма погрешность абсолютная. И вообще, лучше её назвать другой буквой рассчитывать с другим знаком:

$$\widehat{y}=y+\Delta y$$

$$err=\widehat{y}-\ln \left( x-\widehat{y} \right)=y-\ln \left( x-y \right)+\Delta y-\ln \left( 1-\frac{\Delta y}{x-y} \right)\approx \Delta y+\frac{\Delta y}{x-y}\approx \Delta y$$

А на комплексной плоскости, как я уже писал, функция многозначная. Тут надо определиться, либо как рассчитывать значения на различных листах, либо как её разрезать, чтобы сделать однолистной. Как вообще, поступают в таких случаях?