2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки


Правила форума


Посмотреть правила форума



Начать новую тему Ответить на тему
 
 Обратная функция
Сообщение22.01.2019, 18:27 
Аватара пользователя


26/05/12
1694
приходит весна?
Заинтересовала меня тут одна функция. Она является решением вот этого уравнения: $$\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-функция Ламберта-Эйлера, и не стоит плодить сущности без необходимости, но всё же. Может быть и разложение в ряд Тейлора в общем виде есть? Как вообще рассчитывают значения специальной функции в произвольной точке комплексной плоскости? Как строго доказать положение точек ветвления?

 Профиль  
                  
 
 Re: Обратная функция
Сообщение22.01.2019, 19:56 
Аватара пользователя


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

 Профиль  
                  
 
 Re: Обратная функция
Сообщение22.01.2019, 21:40 
Заслуженный участник
Аватара пользователя


22/06/12
2129
/dev/zero
Если убрать строгость, то рецепт такой. Для функции $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$ не знаю, может быть, вы увидите.

 Профиль  
                  
 
 Re: Обратная функция
Сообщение22.01.2019, 22:05 
Аватара пользователя


26/05/12
1694
приходит весна?
А верно, что точки ветвления обратной функции связаны с точками, в которых производная прямой функции обращается в ноль? Например, синус и косинус равны $\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$, то есть как раз те, что у меня получаются из численного счёта.

 Профиль  
                  
 
 Re: Обратная функция
Сообщение22.01.2019, 22:32 
Заслуженный участник
Аватара пользователя


22/06/12
2129
/dev/zero
B@R5uk в сообщении #1370887 писал(а):
вот этого уравнения: $$\exp\left(y\right)+y=x$$

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

-- 22.01.2019 в 22:36 --

StaticZero в сообщении #1370954 писал(а):
для $x \geqslant 0$ не знаю, может быть, вы увидите.

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

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

 Профиль  
                  
 
 Re: Обратная функция
Сообщение23.01.2019, 02:03 
Заслуженный участник
Аватара пользователя


22/06/12
2129
/dev/zero
Процесс
$$
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}$.

 Профиль  
                  
 
 Re: Обратная функция
Сообщение25.01.2019, 13:55 
Аватара пользователя


26/05/12
1694
приходит весна?
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}$) за три итерации. Можно ли сделать это для всей комплексной плоскости?

 Профиль  
                  
 
 Re: Обратная функция
Сообщение25.01.2019, 14:31 
Заслуженный участник
Аватара пользователя


22/06/12
2129
/dev/zero
B@R5uk в сообщении #1371713 писал(а):
То есть лучшего и желать не надо.

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

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

 Профиль  
                  
 
 Re: Обратная функция
Сообщение25.01.2019, 17:19 
Аватара пользователя


26/05/12
1694
приходит весна?
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$$

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

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 9 ] 

Модераторы: Модераторы Математики, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group