Подходы к численному решению уравнения

[Бодигрим]

Подскажите, пожалуйста, с какой стороны подступиться к построению численного решения уравнения

$$x={\int_{-1}^1 p(t)e^{{a\over\tau}xp(t)}\,dt \over \int_{-1}^1 e^{{a\over\tau}xp(t)}\,dt},$$

где $p(t)=(3t^2-1)/2$, $a$ и $\tau$ - некоторые константы. Я не нашел ничего похожего в своих учебниках. Разве что выбрать начальное $x$ и применять к нему последовательно преобразование из правой части - но я слабо представляю себе, как оценить сходимость такого процесса.

[Алексей К.]

Мне в вашем вопросе увиделась относительно простая задача: решить численно уравнение

$$x=\dfrac{f(x)}{g(x)},$$

где $f(x), g(x)$ --- вроде как вполне численно считаемые функции.
Даже всякие умные Maple-Labы захотят в это вопросе пределов изменения $x$. Т.е. будут просить, но неназойливо.

Мне сейчас трудно предположить-поисследовать (потому что ночь, завтра будет трудно, потому что уже не выходной), может можно что-то представить о поведении ${f(x)$ и $g(x)}$. Но, полагаю, не мне одному для дачи более толкового ответа нужно больше подробностей (на форумах Lingvo говорят типа "что ты одну фразу привёл, контекст давай!") --- какие трудности?
Полагаю, реализация самого тупого решения --- "строим график правой части" --- у Вас затруднений не вызывает? Зависимость правой части от $x$ (а вдруг монотонна!) поисследовалась? Ничо не вышло?

[Gafield]

По существу, здесь один параметр, заменим $a/\tau$ на $a$.
$x=0$ - корень для любого а. При $a<4.487657...$ он единственный. При $a=4.487657...$ появляется двойной корень $x\approx0.3236$. При $a>4.487657...$ три корня. Один из из них смещается вправо с ростом $a$, а другой влево и при $a=5$ проходит через ноль, в этот момент корней становится два. А затем снова три. При $a\to+\infty$ левый корень стремится к $-1/2$, а правый к $1$. И все это проще всего заметить, как уже сказано Алексеем К., действительно строя графики разности левой и правой частей. А заметив, уже доказывать. Например, интегралы в числителе и знаменателе выражаются через ierf, поэтому раскладывая в ряд Тейлора правую часть в нуле по х, можно убедиться, что при значении параметра $a$ равном в точности 5 корень $x=0$ будет двойным. Предельные значения корней тоже имеют простое объяснение.

[Бодигрим]

Ага, что-то начало вырисовываться. Спасибо!

Я свел исходное уравнение к степенному ряду с коэффициентами $\asymp 1/k!$. Теперь естественно было бы перейти от ряда к конечной сумме и разрешить полученное алгебраическое уравнение. А как определить возникающую при обрезании ряда погрешность, если границы на $x$ уже определены? Есть какие-то классические результаты?

[Gafield]

Зачем так мучаться? Используя сказанное выше, нахождение корней для конкретных значений $a$ очень просто. Если $a\in(-\infty,a_0)$, где $a_0=4.487657...$, то считать нечего, единственный корень $x_1=0$. Если $a>a_0$, то правый корень $x_2$ ищется на отрезке [0.3235,1], а левый $x_3$ - на отрезке [0,0.3236] или на отрезке [-0.5,0] в зависимости от того меньше ли $a$ чем 5 или больше. А для нахождения единственного корня на отрезке в математических программах типа Maple или Mathematica есть стандартные процедуры. Ничего не надо вычислять вручную. Разве что когда $a$ близко к $a_0$, то прижется уточнить $a_0$.

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

[Бодигрим]

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

[Gafield]

На чем программировать? Разве нет в современных средах типа C++ или Delphi функции интеграла ошибок erf? А erfi - ее разновидность. Вот как у меня посчиталась правая часть в Mathematica:

$$ \frac{1}{2} \left(-1+\frac{e^{\frac{3 a x}{2}} \sqrt{\frac{6}{\pi }}}{\sqrt{a} \text{erfi}\left(\sqrt{\frac{3}{2}} \sqrt{a} \sqrt{x}\right) \sqrt{x}}-\frac{1}{a x}\right) $$

Если нет среди стандартных функций, вероятно, есть в каких-нибудь библиотеках.
Или это учебная задача?

[Бодигрим]

Да, учебная. Иначе не было бы никаких проблем.

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

[Gafield]

Решать через алгебраические уравнения сложное дело. Спрашивается, если решается уравнение с погрешностью $\varepsilon$: $p(x)=\varepsilon$ вместо $p(x)=0$, то с какой точностью находятся его корни? По моему, лучше использовать метод деления отрезка пополам.

Для многочлена Тейлора ошибка оценивается через максимум следующей производной на отрезке на некоторыю константу.

Кстати, а какая функция раскладывается? Не правая ли часть в нуле? Уравнение имеет вид $x-H(ax)=0$. Его можно заменить на уравнение $y/a-H(y)=0$ и в разложении правой части параметра не будет. Первые коэффициенты:

$$ H(y)=\frac{y}{5}+\frac{y^2}{35}-\frac{y^3}{175}-\frac{y^4}{385}+\frac{23 y^5}{875875}+O\left(y^6\right). $$

Экспериментально, первые сто коэффициентов будут убывать степенным образом, что- то вроде $0.26^{n}$. Если такая асимптотика верна, получается, что ряд сходится где-то при $|x|<4$ и использовать его для нахождения корней (по крайней мере, для всех $a$) не получится.

[Бодигрим]

Цитата:

Для многочлена Тейлора ошибка оценивается через максимум следующей производной на отрезке на некоторыю константу.

Нет, я немного не о том. Пусть

$$\sum_{k=0}^\infty {c_k x_0^k\over k!}=0$$

и

$$\sum_{k=0}^n {c_k x_1^k\over k!}=0$$

. Что можно сказать про $|x_0-x_1|$?

Цитата:

Если такая асимптотика верна, получается, что ряд сходится где-то при $|x|<4$ и использовать его для нахождения корней (по крайней мере, для всех $a$) не получится.

Может и впишемся, на самом деле в исходных обозначениях $a/\tau \in [3.78;9.08]$.

У меня, кстати, коэффициенты вышли более быстро падающими. Надо будет завтра перепроверить выкладки...

Спасибо за помощь!

[Gafield]

В общем случае ничего. Например, при $f(x)=\cos(x)$ и достаточно большом $n$ многочлен Тейлора $T_n(x)$ имеет корень $x_1\approx\pi/2$, а $x_0=\pi (k+1/2)$, $k\in\mathbb Z$.

Если предположить, что точки $x_0$ и $x_1$ близки, то можно получнить оценку при некоторых предположениях. Обозначим $r_n(x)=\sum_{k=n+1}^\infty c_k x^k/k!$ остаточный член, $f(x)=T_n(x)+r_n(x)$, тогда, если $r_n(x_0)\ne0$, по теореме о среднем

$$ |x_1-x_0|=\left|\frac{r_n(x_0)}{T'_n(\theta)}\right|, $$

где $\theta$ - некоторое число между $x_0$ и $x_1$. Если в той окрестности $x_0$, где ищется корень, $|f'(x)|>s>0$, и $|r'_n(x)|<s/2$, то получим

$$ |x_1-x_0|\le\frac{|r_n(x_0)|}{|f'(\theta)|-|r'_n(\theta)|}\le2\frac{|r_n(x_0)|}{s}\right|. $$

Если корень простой, $f'(x_0)\ne0, $такую окрестность подобрать можно. Только при $a_0$ корень будет двойным, а если $a$ близко $a_0$, то производная в корне будет мала.

Не понятно, чем метод сведения к алгебраическому уравнению будет в данном случае лучше деления отрезка пополам. Уравнение еще решать надо, и, например, если $a$ чуть меньше $a_0$, то у исходного уравнения нет корней, кроме нуля, а у $T_n$ не исключено, что могут тоявиться какие-то лишние корни. Разве что, посчитать $a_0$ с большой точностью, чтобы отсекать этот случай.

[Бодигрим]

Спасибо за разъяснения про многочлен Тейлора.

Цитата:

Не понятно, чем метод сведения к алгебраическому уравнению будет в данном случае лучше деления отрезка пополам.

А уравнение потом все равно придется дихотомией решать. Но вычисление многочлена куда проще, чем вычисление интеграла.

[Gafield]

Можно преобразовать выражение $x/a-H(x)$, перенеся erfi в числитель и искать нули такой функции:

$$ \frac{\sqrt{\frac{\pi }{6}} \text{erfi}\left(\sqrt{\frac{3}{2}} \sqrt{x}\right)}{\sqrt{x}}-\frac{a e^{3 x/2} \sqrt{\frac{6}{\pi }}}{2 x^2+a x+a}. $$

Первое слагаемое разлагается в ряд с коэффициентами, которые вычисляются явно (почленным интегрированием экспоненты) и убывающими как факториал. Только с выбором интервалов надо быть аккуратнее, при $x<0$ возможен полюс.

[Бодигрим]

Да, я примерно так и сделал в своем решении.