Обобщенная функция Миттаг-Леффлера

Где-то года 2 назад в RTFM_Helpers, я просил помочь вычислить сабж, но...

Функция имеет вид

Код (Text):

1.  
2.                          an
3.        a    (8)     a   x
4. E   (-x ) = SUM (-1) -------
5.  a,b        n=0      Г(ax+b)
6.  

где Г(z) = (z-1)! для z целых, гамма-функция Эйлера.

Так вот при х(0, ~12*pi), все в порядке, а при x>12*pi фунция не вычисляется, точнее вычисляется, но не правильно.

Функция облалает многими интересными свойствами, и при некоторых a,b она превращается в обычные напр. exp(x), sin(x), cos(x), и т.п.

Вычисления проводил обычным суммированием, следует полагать, что необходимо придумать какой-то другой способ более умный. Это и вопрос.

 

http://www.mccme.ru/free-books/lvo/can.html

гляньте сюда, может поможет?



вообще, мысль такая: может этот ряд плохо сходится при указанных значениях х. Если бы ряд сильно расходился, то можно было бы выкрутиться - расходящийся ряд часто ведет себя так - сначала быстро сходится, а потом улетает длеко и надолго, т.е. начинает расходиться.



как считать - честно сказать, не знаю. даже не слышал. сам занимаюсь другим, но поройтесь в Сети - может есть где.

 

Сабж абсолютно сходится при любых х при а!=0.

Проблема в сложении очень больших и очень маленьких чисел как-то надо избавится от этого.

 

Есть мысль перейти к произведению :

Код (Text):

1.  
2.                                   an  
3.                              n   x
4.         a        / (8)   (-1) -------  \
5. E   ( -x ) = Log | MUL 2       Г(an+b) |
6.  a,b            2\ n=0                 /
7.  

надо проверить.

 

murtix



Может быть имеет смысл использовать LaTeX-овскую запись, а то очень уж сложно читать формулы в таком виде. И прочитав нет гарантии что понял правильно..

 

LaTeX -ом пользоваться я к сожалению не умею пользоваться, но могу в таком виде:



565437819__Безымянный.GIF

 

В прошлом году в Stuttgarte был семинар на похожую тему. Посмотрел сейчас на страницу того товарища, который его читал:

http://www.ica1.uni-stuttgart.de/~hansjoerg/

Там сказано(в самом низу) что его работа Numerical Results for the Generalized Mittag-Leffler Function принята к публикации. Так что можно посмотреть этот журнал, вдруг уже напечатали. Или если очень нужно связаться напрямую с автором.



Если есть желание разобраться своими силами(то есть силами форума ) то пока не хватает информации. Сходу возникает довольно много вопросов: как вычислял значение гаммы-функции, для каких значений параметров неправильно вычисляется, откуда "правильные" значения для проверки etc. Опиши подробно что и как вычислял(лучше всего конечно код), попробую воспроизвести проблему у себя и если получится будем думать дальше, но уже предметно.



slow

В принципе было бы действительно неплохо оценить сходимость(т.е. её скорость), чтобы не получилось как с рядом 1/n Но если до 12pi все было нормально, то должно работать и дальше, судя по функции вряд ли 12pi здесь какая-то критическая точка, хотя конечно могу и ошибаться. И опять же не известны параметры.



murtix

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

 

Спасибо за участие.

В общем предметно собственно программа, код и вырезка из моих формул в (*.doc). Только на GDI+, если не ХР, то "gdiplus.dll" могу выслать на мыло, ну чтоб не искали.

Запустите Calculate->Mittag-Leffler и в диалоге просто жмем ОК (диалоговы заморочки потом доделаю). Видим график: в конце он уходит вниз и это неправильно.

Гамма-функцию вычисляю по апрокс. формуле точность отличная проверял факториалом, чем больше аргумент тем точнее.

Насчет сходимости: Гамма-функция сильнее любой степенной функции так что ряд сходится абсолютно.

Насчет "правильных" значений дело в том что эта функция при некоторых a, b переходит в элементарные, например sin(x), его значения-то нам хорошо известны.



5519811__Fm_arh.rar

 

Кому влом запускать прогу и/или морочится с гди+ битмап графика

_798714009__fm.GIF

 

Можно разложить в произведение Г(an+b), и x^an тоже. Потом сгруппировать множители, каждый из них будет не очень большой, я думаю. В сумме получаются слогаемые типа:

х^(an)/Г(an+b) = [x/(an+b-1)]*[x/(an+b-2)]*...*[x/(an+b-k)] * [x^(an-k)/Г(an+b-k)], где an-k<=1,тогда все множители находятся в пределах разумной величины(нужно еще умножить на (-1)^n). Насколько я понимаю, проблема в вычислении не суммы, а слогаемых.

 

Непонятно кстати и почему Е берется от такого странного аргумента: -х^а. Не проще ли произвести замену х:=-х^а и вычислять Е(х)?



С кодом еще не разбирался, сейчас пытаюсь убедить Matlab 6.0 (R12) работать на Pentium 4

 

Нарисовал ту же функцию Е_1.5,1(x) в Matlabe (вычислял суммированием первых 150 членов), кое-что прояснилось. В аттаче кусок для х в интервале [-280,-100], как видно самое интересное начинается как раз после ухода вниз



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



murtix

Таким образом следующий вопрос: откуда известно что это неправильно?



Чтобы проверить алгоритм, вычислил Е_2,1(х) в промежутке от -900 до 0. Результат точно совпал с cos(sqrt(|x|)), как ему и полагается. То есть способ вычисления тут скорее всего ни причем.

 

аттач

1456534114__e1_5_2.jpg

 

>blood

Проблема вычисления именно суммы, а не слагаемых.

Г(х) вычисляется отлично и с х^a тоже никаких проблем, а при умножении и делении не происходит потери точности,

т.е. значащих цифр. А при сложении и вычитании происходит

я в дебуггере видел как 1-1=4,9е-8, понимаешь в чем фишка

Слагаемые данного ряда колеблются

в пределах от 1е+1000 до 1е-18, представь себе всю хрень которая происходит при сложении (и/или вычитании) этих чисел . Именно по этой причине я предлагал перейти к произведению (но тут, как я теперь понимаю, имеются свои заморочки)

>Stiver

Никаких проблем можно перейти к x = -x^a

Просто изначально формула имела такой вид. Если ты заметил в коде я перехожу к x = x^a, я почему-то все время полагал что учитывать знакочередуемость ряда во время вычисления лучше, возможно я ошибался. Посмотрю обязательно.

>Stiver

откуда известно что это неправильно?

Ты же сам ответил на этот вопрос

 

murtix

Честно говоря не заметил.. А можно тогда ответ еще раз?

 

Мда.. дальше в лес, больше дров. Упростить алгоритм вроде бы удалось, т.к. в случае alpha=1.5 мы имеем дело только с "хорошими" аргументами типа n и n+1/2, то можно обойтись без вычисления Г(х). Заодно избавимся и от проблемы больших чисел. (Matlab, так как компактней)

Код (Text):

1.  
2. function res=E_func(x1,x2) % alpha=1.5, beta=1
3.  
4. res=[];
5. for x=x1:0.1:x2
6.     y=1+x/gamma(1.5+1);
7.  
8.     inn=innerk(x,100);
9.     for k=98:-2:2
10.         delta=(x^2)*8/((3*k+3)*(3*k+5)*(3*k+7));
11.         inn=innerk(x,k)+inn*delta;
12.     end
13.     inn=inn*(x^2)*(2^7)/(4*5*6*7);
14.     y=y+inn;
15.    
16.     res=[res [x;y]];
17. end
18. % end E_func ------------------------------------
19.  
20. function r=innerk(z,k)
21.  
22. x=1.5*k+1;
23. r=x/2;
24.  
25. i=1;
26. n=x-1;
27. while 1
28.     r=r*((x+i)/(4*i));
29.     if i==n break    
30.     r=r*((x+n)/(4*n));    
31.     i=i+1;
32.     n=n-1;
33.     if i>n break
34. end
35.  
36. r=r+z/((x+0.5)*sqrt(pi));
37.  
38. % end innerk ------------------------------------
39.  

Используются первые 102 слагаемых, самые большие числа: в процедуре innerk r может быть до 10¹⁰, но и то только в процессе умножения, результат innerk будет всегда около 10. От -240 до 0 функция вычисляется "правильно", то есть совпадает с результатами простого суммирования. Дальше начинаются расхождения. Приведенному выше варианту соответствует картинка e1_5_iter.jpg, видно что график уходит все дальше вниз и поворачивать наверх не собирается. От осциллирования избавиться тоже не удалось.



Интересно, что если в процедуре innerk следующим образом частично поменять i и n местами

Код (Text):

1.  
2.     r=r*((x+i)/(4*n));
3.     ...
4.     r=r*((x+n)/(4*i));    
5.  

то получится картинка e1_5_iter2.jpg, хотя возвращаемое innerk значение естесственно не изменится. Тут я даже не знаю что и думать



В общем пока имеем три разных варианта поведения функции после -240, какой из них правильный непонятно.

1250339213__e_iter.zip

 

>Stiver

Ну при некоторых a и b данный ряд переходит в некоторые элементарные например sin(x)/x, cos(sqrt(|x|)), exp(-x) т.е. в их разложения в ряд, значения этих функций нам известны, но при вычислении обнаруживаем, что функция ведет себя должным образом только на отр. [0, ~12*pi], а дальше фигня. В свое время я изучал поведение членов данного ряда и я уверен что никакой дрожи и тем более ухода вниз быть не должно.

Я даже пробовал способ выч. многочленов

a0+x(a1+x(a2+x(a3+x(...x(an))))), не помогло .

E(a, b, -x^a)->E(1, 1, -x)=exp(-x)

Проверь моей программой

Я дописал "диалоговы заморочки" след. релиз.



_1833805269__Fm_arh.rar

 

>Stiver

Упростить алгоритм вроде бы удалось, т.к. в случае alpha=1.5

Мне жаль тебя разочаровывать, но alpha = alpha(x).

Т.е. Альфа является функцией от х.

 

На самом деле это (alpha = alpha(x)) конечно совсем необязательно, просто я имел ввиду то, что такое возможно.

 

Посмотри пожалуйста на странице

http://mathworld.wolfram.com/Mittag-LefflerFunction.html

вторую сверху формулу, т.е. переход к generalized hyperbolic functions. Или я чего-то не понимаю, или там должно быть ...=Е_a(z^a)(не отсюда ли твой -x^a?) и ограничение только для целых а. Не говоря уже о двух а сразу после "functions", наверное опечатка?



murtix

Поправил: понятно, я просто для примера 1.5 взял, чтобы с твоим графиком сравнить.

12pi здесь ни причем, у каждой фунции этот отрезок свой. У меня Е(1,1,х) пробивает вверх на х=-34, Е(1.5,1,х) на х=-260, а Е(2,1,х) на х=-1200. Чем больше функция похожа на exp(x) тем нестабильнее она вычисляется и быстрее уходит вниз или вверх.