Russian Qt Forum

Товарищи, есть функция вида x = y * exp(y). Как вывести y??  :D

Вывести куда? На улицу и там разобраться? :) Наверное имелось ввиду "выразить y через x", т.е. найти обратную ф-цию. Точно помню - производная обратной равна 1 / производную прямой

(y * exp(y))' = exp(y) * (1 + y)  // производная прямой ф-ции
Значит обратная будет =
Ну конечно такой интеграл я не возьму - но погуглите, помню видел взятие интегралов в online. Не исключено что он берется

Исключено (в элементарных функциях)

намек на то, что выразить y через х не удастся?

В элементарных функциях - нет, не удастся. Доказано. Можно через специальную функцию "интегральный синус"
https://ru.wikipedia.org/wiki/%D0%98%D0%BD%D1%82%D0%B5%D0%B3%D1%80%D0%B0%D0%BB%D1%8C%D0%BD%D1%8B%D0%B9_%D1%81%D0%B8%D0%BD%D1%83%D1%81

Аналитически вы не выразите y через x..

Но численно это делается достаточно просто..

Ну, численно это понятно, график построить всегда можно.
Но вот интегральный синус меня печалит..  :'(

Или намек не просто на график, а на вывод зависимости по графику?

ОКТА, оставьте напрасные надежды. Интеграл exp(x)dx/x в элементарных функциях НЕ БЕРЕТСЯ. Как и все, к нему сводящиеся. Доказано, ИМХО, еще до Эйлера.

Если верить маткаду, то в аттаче решение.
Ei - интегрально показательная функция (https://ru.wikipedia.org/wiki/%D0%98%D0%BD%D1%82%D0%B5%D0%B3%D1%80%D0%B0%D0%BB%D1%8C%D0%BD%D0%B0%D1%8F_%D0%BF%D0%BE%D0%BA%D0%B0%D0%B7%D0%B0%D1%82%D0%B5%D0%BB%D1%8C%D0%BD%D0%B0%D1%8F_%D1%84%D1%83%D0%BD%D0%BA%D1%86%D0%B8%D1%8F)

Да, тут на работе уже кандидат один и доктор подтвердили, что никак  ;D
Но все равно, спасибо всем!  :D

Сейчас вам в уши тут надуют)
То что там выше понаписано - полная ерунда.. я про это:
Вообще непонятно что это даст..

Можно свести это выражение к интегральной форме, но в итоге это всё равно приведёт к начальной неявной зависимости y(x)..

Только вот это решение не имеет никакого отношения к первоночальной задаче..

Да, неверно. Интегрировать-то надо по dx (в не по dy). Т.е. когда получили производную прямой - надо выразить ее через x. Пример

x = exp(y) // прямая
x' = exp(y)  // производная прямой
y' = 1/x' = 1/exp(x) => 1 / x  // выражаем y через x
интеграл(1/x) dx = ln(x)  // нашли обратную

Но в случае y * exp(y) избавиться от y не удается

Если нетрудно, покажите доказуху  :)

Доказухи не нашел. Но вот такая парочка ссылок. Я думаю, что если бы интегральная экспонента выражалась бы через элементарные, то не было бы смысла вводить для нее СПЕЦИАЛЬНУЮ функцию. Но это, конечно, не доказательство, а так, общие соображения
https://ru.wikipedia.org/wiki/%D0%98%D0%BD%D1%82%D0%B5%D0%B3%D1%80%D0%B0%D0%BB%D1%8C%D0%BD%D0%B0%D1%8F_%D0%BF%D0%BE%D0%BA%D0%B0%D0%B7%D0%B0%D1%82%D0%B5%D0%BB%D1%8C%D0%BD%D0%B0%D1%8F_%D1%84%D1%83%D0%BD%D0%BA%D1%86%D0%B8%D1%8F
https://ru.wikipedia.org/wiki/%D0%A1%D0%BF%D0%B5%D1%86%D0%B8%D0%B0%D0%BB%D1%8C%D0%BD%D1%8B%D0%B5_%D1%84%D1%83%D0%BD%D0%BA%D1%86%D0%B8%D0%B8

Вообще какая то раздутая проблема..

Это решается банальным нахождением корня уравнения x = y*exp(y).. Например, обычным bisection..

C++ (Qt)
#include <iostream>
#include <fstream>
 
#include <boost/math/tools/roots.hpp>
 
template <class T>
T my_func(const T & x, const T & err = 0.000001)
{
    using namespace boost::math::tools;
    const T min = 0;
    const T max = x; // x >= 0
 
    auto res = bisect([&](const T& s) { return x - s * exp(s); }, min, max,
    [&](const T& a, const T& b) { return std::fabs(b-a) <= err; });
 
    return (res.first + res.second)/2;
}
 
int main()
{
    double x = 0;
    const double dx = 0.01;
 
    std::ofstream out("data.txt");
 
    while (x < 100)
    {
        out << x << " " << my_func(x) << std::endl;
        x += dx;
    }
 
    return 0;
}
 

При x < 0, решения не существует вовсе.. А вот при каких x решения не существует - вопрос интересней)

А что даст корень уравнения?

Корень и есть "какое y соответствует данному x". Но вот считать его для каждой точки - подозреваю будет накладно. Кстати ф-ция немонотонна и неизвестно правильно ли это вообще.

Да, и возможна ситуация когда одному x соответствует неск y, т.е. результат не определен. Поэтому по-простому-народному выберите диапазон y (напр от 0 до 10), там все монотонно, создайте массив напр 1000 значений x и спокойно находите через lower_bound

Функция как раз монотонна, и однозначна (во всяком случае для x >= 0). Бисекшин - достаточно быстрый метод, никаких "накладно" вы там не заметите..   

А вот так лучше как раз не делать)

Ну что Вы мне парите :'( Пусть бисекшин - простая дихотомия, т.е. деление отрезка пополам (иначе еще хуже). Напр 10 шагов - уже 10 раз вычислена exp, а это разложение в ряд и.т.п.

Аргументы?

Посижу, почитаю, поучусь у умных дядек  ;D

Вы такую задачку не слышали: Если сложить газету 40 раз, то толщина стопки получится равной расстоянию от земли до луны..
Так вот к чему это я? А теперь прикиньте, какую точность мы получим  за всего каких то 10 итераций бисекшином..

А они для Вас не очевидны? Давайте сравним две реализации.. С условием, что вычислять функцию нужно с заданной точностью при любом положительном x. Думаю тогда все вопросы сразу отпадут)
 

Рано меня ещё в "дядьки" записывать)

bisect требует диапазона аргумента (в данном случае y). Откуда Вы его возьмете? Ну напр y = 50, exp(50) - уже "слишком много" чтобы реально использоваться в задаче. Поэтому разговоры о "любом положительном x" неуместны. Более того, на практике часто вообще интересен диапазон [0..1]. Так что давайте реалистичнее, тогда пожалуйста - сравним  :)

Верхнюю границу y можно оценить, и поверьте, она  50ти никогда не будет равна, для любого разумного x..
Так что разговор уместен)

Это Вы с потолка опять взяли? А мне нужно график, например от 0 до 1000000.. И что мне с вашим от 0 до 1 делать? Ммм?

 

Я говорю о значениях аргумента (y). exp(100000..) вернет nan. Наверное Вы имеете ввиду график х(y) где макс x - достаточно большое. Ну что ж, потрачу больше памяти чтобы добиться нужной точности. Да, если Вы не поняли (судя по постам) то массив строится 1 раз на этапе предрасчета. Поэтому соревноваться с ним в скорости поиска довольно глупо  :)

Я понимаю, о чём Вы говорите.. Я тоже говорю о значениях y..  Ситуации с nan никогда не возникнет, поскольку верхняя граница y для всех разумных x, будет порядка 10, а то и меньше.. Я имею в виду функцию y(x)

Давайте, проверим) Выкладывайте свою реализацию "по-народному", а я выложу свою) 
Условие прежние посчитать y(x) для любого положительного x с заданной точностью)

  

Похоже, народное решение мы не увидим.. Ну да ладно)

Скажу несколько слов о производительности метода bisect..
Написал пробный тест, где сравнивается скорость построения графика нашей функции (exact) y(x) (обратной к y*exp(y))
и очень на неё похожей функции (approx) y = log(1 + x)/1.35. На приаттаченом рисунке они обе изображены.

Переменная x проходит от 0 до 10000 с шагом dx = 1. Результат записывается в файл.

Исходник с замером времени:

C++ (Qt)
#include <iostream>
#include <fstream>
#include <chrono>
 
#include <specmath/roots.h>
 
 
template <class T>
T exact_func(const T & x, const T & err = 0.000001)
{
    static const T e = exp(1);
 
    const T ln_x = log(x);
    const T min = (x < e) ? 0 : 1;
    const T max = (x < e) ? 1 : ln_x;
 
    return specmath::bisect([&](const T& s) { return (x < e) ? (x - s * exp(s)) : (ln_x - s - log(s)); },  min, max, err);
}
 
template <class T>
T approx_func(const T & x)
{
    return log(1 + x)/1.35;
}
 
int main()
{
    double x = 0;
    const double xmax = 10000;
    const double dx = 1;
 
 
    std::ofstream out1("data1.txt");
    std::ofstream out2("data2.txt");
 
    auto start = std::chrono::high_resolution_clock::now();
    while (x < xmax)
    {
        out1 << x << " " << exact_func(x) << std::endl;
        x += dx;
    }
    auto stop = std::chrono::high_resolution_clock::now();
 
    auto elapsed = std::chrono::duration_cast<std::chrono::milliseconds>(stop-start).count();
 
    std::cout << " bisection: total time = " << elapsed << "ms" << std::endl;
 
 
 
    x = 0;
    start = std::chrono::high_resolution_clock::now();
    while (x < xmax)
    {
        out2 << x << " " << approx_func2(x) << std::endl;
        x += dx;
    }
    stop = std::chrono::high_resolution_clock::now();
 
    elapsed = std::chrono::duration_cast<std::chrono::milliseconds>(stop-start).count();
 
    std::cout << "aprox func: total time = " << elapsed << "ms" << std::endl;
 
    return 0;
}
 

 

Результаты трёх прогонов следующие:

Bash
 bisection: total time = 58ms
aprox func: total time = 50ms
 
 bisection: total time = 58ms
aprox func: total time = 50ms
 
 bisection: total time = 58ms
aprox func: total time = 51ms
 

Как видно, разница ничтожна..  
"Народный" вариант от igors однозначно проиграл бы по памяти (это мягко говоря), да и по скорости тоже есть сомнения..
Так что делаем выводы..

ЗЫ В примере выше, вместо boost::bisect использовалась более  быстрая её версия.. Но сути это не меняет..

Что-то не то, не так. В approx всего лишь один log, в exact уже до bisect тот же log плюс сам bisect чего-то весит, сколь бы чудесным он ни был. Вычислительная нагрузка в exact по меньшей мере в неск раз больше. Чего же это времена отличаются ненамного? Напечатайте сколько раз зовется функтор bisect, тогда будет ясно

f = ln_x - s - log(s)  // прологарифмированное исходное
f(1) = log(1) - 1 - log(1) = -1  // нижний предел
f(log(x)) = log(x) - x - log(log(x)) < 0 // верхний предел

Похоже bisect'у нечего считать т.к. оба предела < 0

Я найду время, просто "не сию минуту" - есть и другие дела  :)

Ну вообще то, эти результаты: два графика были построены после запуска этого примера.. Строится всё корректно) Время тоже считается честно)  

Нет, если внимательно посмотреть, то в самой лямбде тоже проверяется условие (x < e) и в зависимости от него выбирается своя схема..
Там всё корректно, иначе результат (график) не получился бы)  

Поменял местами exact и approx.. Результаты не изменились..

Приаттачиваю проект, на всякий случай..
 

Да-да-да. Вам же нужно еще показать нам "профессиональную" реализацию (http://www.prog.org.ru/index.php?topic=26959.msg194871#msg194871) - просим просим. :)
Да и неплохо бы прокомментировать свои традиционно голословные заявления по поводу невозможности сделать прокси устройство от QIODevice, которое умеет испускать исключения.
Есть еще "другие дела", есть.

А то я уже начал расстраиваться... думаю забыл мастер сирых. Ан нет.  ::)

P.S. И название темы подобрано хорошо: Мать И Матика. :)

Ага, понял, у Вас все время уходит на запись в поток. В аттаче обещанная реализация. Насчет точности, расхода памяти, удобства - тут и спорить нечего, bisect превосходен. Но все-таки табличное задание быстрее, и может быть еще быстрее если ф-ция сложнее. И этот плюс может перевесить все минусы.

Если Вы в припадке раздражения - то я лучше промолчу, все равно разговор ничем хорошим не кончится  :)

Вы меня расстраиваете, "профессиональную" сортировку можно не ждать? :(
Все как всегда, пришел, бла-бла-бла и в кусты. Зато как щеки раздувал. Дааа, это не финдреплейсы на QStringах писать... понимаю. :)

О чем Вы ??? Какую сортировку я Вам обещал  ???
(классика советского кино)

Ну как же "профессиональную". :)
Ну понятно, ниасилил, бывает.

Посмотрел) Надо отметить - изысканная вещь для ценителей хардкора)
Я не представляю себе пользователя, который бы использовал это в реальном проекте..
Да, в этом конкретном случае, работает быстрее, но учитывая, что для  реальных задач важна не только скорость (а скорость в обеих случаях неплохая и едва ли будет где-либо узким местом) но и точность (и память) то оно того не стоит..
Таким решением просто неудобно пользоваться, к тому же оно ещё и мой бисекшн использует :)

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

Поэтому такое решение имеет очень ограниченный круг применения..

Кстатии, эту функцию (обратную к y*exp(y)) можно найти приближённо аналитически (по теория возмущений):

C++ (Qt)
template <class T>
T approx_func2(const T & x)
{
    static const T r = 1.33;
    const T y0 = log(1 + x)/r;
 
    const T a = 0.5 / (y0*y0);
    const T b = -(1 + 1.0/y0);
    const T c = log(x/y0) - y0;
 
    return y0 + (-b - sqrt(b*b - 4*a*c))/(2*a);
}
 

С абсолютной погрешностью меньше 0.0003 на всём интервале x от 0 до нескольких миллионов..

В общем,  мораль, как всегда, такова: лучше вначале хорошенько 10 раз подумать (подальше от компьютера) и только потом уже кодить  ;) 


 

 
   

Это для школьников, которые лабы делают.
А настоящие "профессионалы" во время звонка клиента, который говорит "Слушай, хочу такую прогу....", уже набирают в редакторе "class MainWindow". Это настоящие мужики. :)

Я бы сказал, это очень брутальные мужики :)

/* Чёта почувствовал себя опять школьником) */

Насколько? Я не знаю что Вашей машине, на моей в 8-9 раз быстрее

Все эти "стоны придавленного дерева" никто не будет даже слушать - если сказано хотя бы "вдвое быстрее", а о большем неудобно и говорить. В задачах рендера за это продают душу (ну у кого она есть)

Вот Вы считаете что я тупой ретроград  :) Ну ладно, пусть так (на Вас не обижаюсь). Но тогда поясните почему я (такой тупой) применил куда более общее решение чем Вы, такой современный и продвинутый? :) Ну в самом деле - то Вы прологарифмировали уравнение, то Вы заменили его приближенным. Это прекрасно - но только в частном случае. А в табличном решении пофиг какая ф-ция - все равно это предрасчет выполняемый один раз. Имеет др минусы? Не спорю, но как и всякое решение (в том числе и Ваше). Кстати в игровых задачах фаза предрасчета может длиться неск суток - и это норм (все равно ноль в runtime)

Не надо судить "по одежке" (типа сколько там template  :)) - это очень поверхностный взгляд на вещи

А Вам вообще "не дано" придумать свое решение - слишком долго сидели на дусте. Ну что ж, "прекрасный реализатор" - тоже неплохо, кто знает...

Ну вообще пока, только вы ничего кроме финдреплейса сделать не осилили. :)

Как в том анекдоте
Может я и был неправ (не судите - и не судимы будете), но в общем хорошо - пока так пока.

У меня примерно также.. Но..

Мне частенько приходиться писать числодробилки, в которых такие трансцендентные уравнения возникают постоянно, и функции там потяжелее логарифмов с экспонентами) Но вот табулировать их значения - это дохлый номер.. И никто (кроме, наверное Вас) этим заниматься не будет по следующим причинам:

Во-первых, по мимо аргумента функции, имеются ещё и доп параметры.. Даже если взять пример TC и чуть усложнить его:
x = y*exp(a * y), то Ваш табличный метод уже в пролёте.

Во-вторых, в данном примере, Вам просто повезло с функцией: y - имеет очень медленную логарифмическую зависимость и поэтому для очень большого интервала x, y меняется незначительно..
А если взять хотя бы такую функцию: x = pow(y, 1/4) + pow(y, 1/6), то размер Вашей таблички очень быстро вырастет если мы претендуем хоть на какую то точность, пусть хоть и во втором знаке..

Нет, это Вы преувеличиваете)

Так о том и речь, что Ваше решение очень ограниченное) Неужели это до сих пор не очевидно?   
 

А почему сами судите?
Вы не стесняетесь судить в темах 8Observer8, не замечали? А в другой теме, вы назвали комментатора мудаком и посоветовали не обращать на него внимание. Забыли?
А когда вас просят показать себя, вы сдуваетесь и начинаете истерить и говорить "сам дурак".
Традиционно, чем меньше вы в теме разбираетесь, тем сильнее вам хочется её комментировать, а когда вам говорят о вашем не знании - начинаете капризить. Ну сколько можно?

Да какое оно "мое"? :) Табличный метод известен с незапамятных времен, и также известно что он ощутимо быстрее аналитики. Хотели проверить, ну вот проверили.

Насчет ограничений никто не спорит. Напр atan табличкой не сделать (можно через cos дважды). Но аналитическое решение тоже имеет минусы. Вот напр пользователь задал точки сгладив их по Безье - ох и долго Вам придется возюкаться с формулами :) А искать "делением пополам" - так монотонности никто не обещал.

Ну... спасибо всем большое за советы и обсуждение  :)
Сейчас начну пробовать ваши идеи на своей задачке и выскажусь по результатам!  ;)

m_ax, а что это за specmath? Путался нагуглить и ничего не нашел  ???

Ну это чисто моё, хэнд мэйд творчество) Набор тех средств, которыми приходится часто пользоваться.. 

Огонь!) Буду нагло пользоваться!!  :D

А скажи пожалуйста, откуда в linear_interpolate два последних условия берутся и для чего они?

Это для устойчивости.. Дело в том, что если разность fb-fa окажется большой, а интервал b-a очень маленьким, то из-за ограничения машинной точности, в результате можете ерунду получить.. Вот чтобы этого не было мы подстраховываемся теми двумя условиями..)

Спасибо за разъяснения, m_ax!  ;)
Не хочу никого обижать, но способ, который предложил m_ax, в условиях текущей задачи оказался предпочтительнее.

Если на нескольких страницах этой темы ответ уже был - прошу прощения, мог пропустить. А так все в меру элементарно. Решения не существует при x<-1/e :)