УДК 519.688
ОСОБЕННОСТИ ВЫЧИСЛЕНИЯ ФУНКЦИИ ТИПА МИТТАГ-ЛЕФФЛЕРА В СИСТЕМЕ КОМПЬЮТЕРНОЙ МАТЕМАТИКИ «MAPLE» Паровик Р.И.1,2
1 Институт космофизических исследований и распространения радиоволн ДВО РАН, 684034, Камчатский край, с. Паратунка, ул. Мирная, 7
2 Филиал Дальневосточного Федерального государственного университета, 683031, г. Петропавловск-Камчатский, ул. Тушканова, 11/1
E-mail: romanparovik@gmail.com
В работе предложены вычислительные алгоритмы для корректного вычисления функции типа Миттаг-Леффлера с дробным параметром 1 < а < 2 в математическом пакете MAPLE.
Ключевые слова: фрактальный параметрический осциллятор, функция типа Миттаг-Леффлера, контур Ханкеля, интегральное представление
(с) Паровик Р.И., 2012
MSC 65R10
CALCULATION SPECIFIC FUNCTIONS OF MITTAG-LEFFLER IN THE COMPUTER MATHEMATICS «MAPLE» Parovik R.I.1,2
1 Institute of Cosmophysical Researches and Radio Wave Propagation Far-Eastern Branch, Russian Academy of Sciences, 684034, Kamchatskiy Kray, Paratunka, Mirnaya st., 7, Russia
2 Branch of the Far Eastern Federal State University, 683031, Petropavlovsk-Kamchatsky, Tushkanova st., 11/1, Russia
E-mail: romanparovik@gmail.com
In this work the computational algorithms to correctly calculate the functions of Mittag-Leffler with fractional parameter in mathematical package MAPLE.
Key words: fractal parametric oscillator, function of Mittag-Leffler, circuit Hankel, integral representation
(c) Parovik R.I., 2012
Введение
В настоящее время широкое развитие получило математическое моделирование фрактальных процессов. Под фрактальными процессами понимаются различные процессы природных, экономических, историко-социальных явлений с учетом их фрактальных свойств. Математическая интерпретация таких явлений дается, хорошо разработанным за более чем 300 лет, аппаратом дробного исчисления [1].
Особенностью математических построений фрактальных моделей является обобщение дифференциальных уравнений целочисленного порядка на вещественный порядок, что приводит к решениям со степенной зависимостью. В качестве такого решения часто выступает целая специальная функция типа Миттаг-Леффлера.
Вычисление этой функции по определению в системах типа «MAPLE» приводит зачастую к некорректным результатам: либо совсем не вычисляется, либо из-за погрешности вычисленний дает неверный ответ. Поэтому в работе мы остановимся подробно на вычислительном алгоритме этой функции, но сначала приведем некоторые ее свойства.
Функция типа Миттаг-Леффлера и ее некоторые свойства
В работе [2] Г. Митагг-Леффлер ввел в рассмотрении обобщенную показательную функцию:
Здесь Г(х) - гамма-функция Эйлера. Из соотношения (1) видно, что если г = 0, то функция Миттаг-Леффлера становится константой Еа (0) = 1, а также в случае а = 1 становится экспонентой Е1 (г) = ег. Более подробно о свойствах функции Миттаг-Леффлера можно узнать в работе [2].
В работе [3] исследуется функция типа Миттаг-Леффлера:
Функцию (2) принято называть функцией типа Миттаг-Леффлера. Она является решением многих дифференциальных уравнений в производных дробного порядка. Например, фрактальное осцилляционное уравнение [4], а также в работе автора [5]:
где С1 и С2 - константы интегрирования, которые определяются из начальных условий. Отметим некоторые свойства функции типа Миттаг-Леффлера, вытекающие из определения (2):
(1)
(2)
d0tu (т) + юв u (t) = О, 1 < в < 2
имеет решение
u (t) = CiEe,1 ( — (®t)в) + C2tEe ,2 ^ — (®t)в) ,
1) Ea,в (О) = ^, Ea,1 (z) = Ea (z) ^ EM (z) = ez;
2) E1,2 (z) = —, E2,1 (z) = cA (V) и E2,2 (z) =
zz
1 dE™ я+1 (z)
3) Ea ,в (z) = ГЩ + ZEa ,в+а (z) и Ea ,в (z) = вEa ,в+1 (z) + Za-^-;
/ d \ m
4) (^J (zв—1Ea,в (za)) = zв —m—1Ea,в —m (za), m > 1;
1 m—1 (
5) Eaв (z) = - Е £а/mJJ (z1/mei2nk/m
m к=О v
6) /zв lEa,в (Яtа) tв 1dt = z^,в+1 (ЯZa).
О
Также справедливы асимптотические разложения функции типа Миттаг-Леффлера при больших значений |z|:
а) при |arg(z)| < а1п, а1 е (а/2,а), а е (0,2), n е N имеем
1 1 — в / 1 \ П z k / 1
Еа ,в(z) = аz "-~exp (z s) - Ег ™_а.л + O
а V J k=lг (в — аk) Vlz|n+1
б) при а1п < |arg(z)| < п, а1 є (а/2, а), а є (0,2), n є N имеем
Eа,в (z) Е г (я - аk) + O I Uin+1 I ;
в) в случае а > 2, n є N
Еа,в (г) а ехр(гт) Е Г(в-ак) + °( |г|”+>
суммирование для первого слагаемого производится для тех т, при которых выполнено условие:
|а^(г) + 2пт| < е + ап/2, е > 0.
Доказательства асимптотических разложений функций типа Миттаг-Леффлера можно найти в работе [3].
Представление функции типа Миттаг-Леффлера в комплексной плоскости
Будем считать в > 0. Рассмотрим контур (петлю) Ханкеля у(е,8) (см. рис.), который состоит из следующих трех частей: луча ^ = {а^(£) = —8, |£ | > е}, луча S-8 = {а^) = 8, |£ | > е} и дуги окружности С8 (0,е) = {—8 < а^(£) < 8, |£ | = е}.
В свою очередь контур у(е, 8) разбивает плоскость £ на две бесконечные подобласти: ) (е, 8) и 0(+) (е, 8). В работе [3] показано, что если выполняются условия
5З
О < а < 2 и ап < й < min{п,па}, то имеют место следующие интегральные представления:
1 г eE а е г—я
Еа* (z) = 2-п J -Ї—ГdE, z є G<-)(£,й), (З)
Y <є,й)
4 кЬЯ
"'"'''(є,
e-,в (z) =1 z 1-в ez 11 2-Гп7 / ^Т-Тdе ,z є oW(e,й)
Y <є,й)
Интеграл, стоящий в соотношениях (3), вычисляется по каждой из частей контура у(е, 5), считая ^ = еег5 и следуя результатам работ [6], [7]:
ап
Еа,в (z) = JK(а,в, й, г,z)drI J P (а, в, є, ф,z)dф, z є G<—^ (є, й)
ап
1 1-в 1 Г
Еа в (z) = аzezа I K (а, в, й, г, z)drI а
ил
I j P (а, в, є, ф, z^, z є (є, й), (4)
ап
P (а, в, £, Ф,z)d0, z е
ап
к (а, в, в, r,z) = -r1-4^-(в/а)r si" (*- 8) ~z sin (2),
ап r2 — 2rz cos (8) + z2
£ £1+(1-в)/а ££®С0Кф/а) (cos (ю) + isin (ю))
(а,в, £, ф,z) - — £егф - z ,
Ф — ri sin (5/а) + 8 (1 + (1 - в)/а), ю — £« sin (ф/а) + ф (1 + (1 - в Vа).
Интегральное представление (4) состоит из двух частей (интегралов): первый интеграл - монотонная часть, а второй - осцилляционная. Рассмотрим случай 0 < а < 1. Тогда получим, что 8 — min{п,ап} — ап. Поэтому будем иметь:
K (а, в, ап, r, z) — K (а, в, r, z) —
ОО
l-в -
r - e ra rsin (п (l - в)) - zsin (п (l - в + a))
ап
r2 - 2rzcos (ап) + z2
Здесь возможны три случая:
1 ) |эд^(г) | = ап, 2) |аг§(г) | > ап, 3) |эд^(г) | < ап.
В работе [7] доказаны теоремы для этих случаев. Приведем результат.
В первом случае функцию типа Миттаг-Леффлера можно вычислить с помощью следующего соотношения:
ап
Еав (г) = J К (а, в, г, г)^г + ^ Р (а, в, £, Ф, г)#, е > |г|. £
Во втором случае
(5)
ап
/КГ (а, в, r, z)dr, в < l + а
Ea,в (z) = < -
sin (ап)
1 — г а
ап
r2 - 2rzcos (ап)+ z2 z
dr —, в — l + a .
(6)
ап
/К(a,в, r,z)dr + / P (a, в, є, ф,z)dф, є > О, в > О v є —ап
В третьем случае
°° _ z( 1-в )/а ez -
f К (а, в, r,z)dr +---------------------, в < l + а
є а
Еа ,в (z) = < -
sin (ап)
1 -Г а
ап
r2 - 2rzcos (ап) + z2
dr +
a z
l -в -
_ ап z - ez
f К (а, в, r, z)dr + f P (а, в, є, ф, z)d0 +----------------
k. є -ап a
в = l + a
(7)
,z a
є > О, в > О
Алгоритм вычисления функции типа Миттаг-Леффлера
В интегральных представлениях (5), (6), (7) интерес представляет несобственный интеграл. Его можно представить в качестве определенного интеграла. Введем, согласно работе [7], фиксированную константу q = 0,9. Рассмотрим следующие три случая:
1 ) |г| < q, 0 < а, 2) |г| > q, 0 < а < 1 , 3) |г| > q, 1 < а.
В первом случае, согласно (2):
кО
Ea,в (z) = Е
к=О
Г (в + a к)
+ П (z), |П (z)| < P,
ко = тах{ [(1 - ^)/^ + 1, |_1п (р(1 -|г|)) / 1п (|гШ} . Здесь р - точность вычисления, Ц - целая часть числа.
оо
оо
e
к
z
Во втором случае
°° Г0
JK(а,в,r,z)dr = JK(а,в,r,z)dr + n (z), |n (z)| < Р, а е (0,е),
а а
Го = max jl, 2 |z|, (-ln (пр/б))а j .
В третьем случае, с учетом свойства 5, где m = [аJ + 1 и 0< а/т < 1 функция типа Миттаг-Леффлера вычисляется по формуле:
1 т-1 / ч
Еа ,в (z) = т Е Еа/т,в (z1/mei2n‘/m
т k=0
для случаев |z|т < ди |z| т > q.
Вычисление функции типа Миттаг-Леффлера с помощью системы компьютерной математики MAPLE
Положим р = 10-4. с учетом методики рассмотренной выше и свойств функции типа Миттаг-Леффлера составим на языке MAPLE следующую процедуру ML1:
> ML1:= proc (alpha, beta, x) local k0, phi, r, rho, r0; r0:=evalf(max(1, 2*abs(x), trunc(-ln(0.0001*Pi/6))‘alpha));
if x=0 then 1/GAMMA(beta) elif alpha=1 and beta=1 then evalf(exp(x))
elif alpha=1 and beta=2 then evalf((exp(x)-1)/x) elif alpha=2 and beta=2 then evalf(sinh(sqrt(x))/sqrt(x)) elif alpha = 2 and beta = 1 then evalf(cosh(x‘(1/2))) elif 0<alpha and alpha < 1 then if abs(x)<0.9 then
k0:=max(trunc((1-beta)/alpha)+1, trunc(ln(0.0001*(1-abs(x)))/ln(abs(x)))); sum(x‘k/GAMMA(beta+alpha*k), k=0 .. k0) elif abs(x)<trunc(10+5*alpha) then if alpha*evalf(Pi)<abs(evalf(argument(x))) and 0.0001 < abs(evalf(abs(argument(x)))-alpha*evalf(Pi)) then if beta<1+alpha then evalf(Int(r‘((1-beta)/alpha)*exp(-r‘(1/alpha))*(r*sin(Pi*(1-beta))-x*sin(Pi*(1-beta+alpha)))/(alpha*Pi*(r‘2-2*x*r*cos(alpha*Pi)+x‘2)), r=0..r0)) else evalf(Int(r‘((1-beta)/alpha)*exp(-r‘(1/alpha))*(r*sin(Pi*(1-beta))-x*sin(Pi*(1-beta+alpha)))/(alpha*Pi*(r‘2-2*x*r*cos(alpha*Pi)+x‘2)),r=1..r0))+ evalf(Int(0.5*exp(cos(phi/alpha)*
exp(I*(phi*(1+(1-beta)/alpha)+sin(phi/alpha))))/(alpha*Pi*(exp(I*phi)-x)), phi = alpha*evalf(Pi)..-alpha*evalf(Pi))) end if
elif abs(x)<trunc(10+5*alpha) and abs(argument(x)) < alpha*evalf(Pi) and 0.0001 <abs(abs(argument(x))-alpha*evalf(Pi)) then
if beta<1+alpha then evalf(Int(r‘((1-beta)/alpha)*exp(-r‘(1/alpha))*(r*sin(Pi*(1-beta))-x*sin(Pi*(1-beta+alpha)))/(alpha*Pi*(r‘2-2*x*r*cos(alpha*Pi)+x‘2)), r=0..r0))+ x‘((1-beta)/alpha)*exp(x‘(1/alpha))/alpha
else evalf(Int(r‘((1-beta)/alpha)*exp(-r‘(1/alpha))*(r*sin(Pi*(1-beta))-x*sin(Pi*(1-beta+alpha)))/(alpha*Pi*(r‘2-2*x*r*cos(alpha*Pi)+x‘2)), r=0.5*abs(x)..r0)+ simplify(evalf(Int(0.5*(0.5*abs(x))‘(1+(1-beta)/alpha)*exp((0.5*abs(x))‘(1/alpha)* cos(phi/alpha))*(cos(phi*(1+(1-beta)/alpha)+(0.5*abs(x))‘(1/alpha)*sin(phi/alpha))+ I*sin(phi*(1+(1-beta)/alpha)+(0.5*abs(x))‘(1/alpha)*sin(phi/alpha)))/(alpha* Pi*(0.5*abs(x)*exp(I*phi)-x)), phi=alpha*Pi.. -alpha*Pi)))+ x‘((1-beta)/alpha)*exp(x‘(1/alpha))/alpha end if else evalf(Int(r‘((1-beta)/alpha)* exp(-r‘(1/alpha))*(r*sin(Pi*(1-beta))-
x*sin(Pi*(1-beta+alpha)))/(alpha*Pi*(r‘2-2*x*r*cos(alpha*Pi)+x‘2)),
r=0.5*abs(x)+0.5..r0))+
simplify(evalf(Int(0.5*(0.5*abs(x)+0.5)‘(1+(1-beta)/alpha) *exp((0.5*abs(x)+0.5)‘(1/alpha)*cos(phi/alpha))*(cos(phi*(1+(1-beta)/alpha)+ (0.5*abs(x)+0.5)‘(1/alpha)*sin(phi/alpha))+I*sin(phi*(1+(1-beta)/alpha)+ (0.5*abs(x)+0.5)‘(1/alpha)*sin(phi/alpha)))/(alpha*Pi*((0.5*abs(x)+0.5)*exp(I*phi)-x)), phi=alpha*Pi..alpha*Pi)))+ x‘((1-beta)/alpha)*exp(x‘(1/alpha))/alpha end if
else k0:=trunc(-ln(0.0001)/ln(abs(x)));
if abs(evalf(argument(x)))<0.75*alpha*evalf(Pi) then x‘((1-beta)/alpha)* exp(x‘(1/alpha))/alpha-(sum(x‘(-k)/GAMMA(beta-alpha*k), k=0..k0)) else simplify(evalf(-(sum(x‘(-k)/GAMMA(beta-alpha*k), k=0..k0)))) end if end if end if end proc:
В случае, когда l < a < 2 или m = 2, составим процедуру ML2:
> ML2:=proc(a,b,z) if z=0 then 1/GAMMA(b) elif a=2 and b=1 then cosh(sqrt(z)); elif a=2 and b=2 then evalf(sinh(sqrt(z))/sqrt(z)) elif 1<a and a<2 then 0.5*(ML1(a/2,b,z‘0.5)+ML1(a/2,b,-z‘0.5)) fi;
end proc:
Результаты вычислений функции типа Миттаг-Леффлера
Пример 1. Сначала вычислим функцию типа Миттаг-Леффлера согласно определению (2):
> а1:=0.7;
> Ь1:=1.3;
> z:=0.7;
> t:=time();
сії - 0.7
Ы = 1.3
г = 0.7
г— 3391.446
> E:=sum(z‘k/GAMMA(b1+a1*k),k=0..infinity);
> TIME=time()-t;
TRIE = 20.233
Вычислим эту же функцию с помощью процедуры МЫ:
> И:=Нше();
(Г= 3431.959
> ML1(0.7,1.3,0.7);
2.294030341
> TIME=time()-t1;
TRIE - 0.
Результаты этого примера показывают, что методика, рассмотренная в работе, позволяет вычислять функцию типа Миттаг-Леффлера с сингулярностями.
Пример 2. Пусть 1 < а < 2.
> а1:=1.3;
а! ■= 1.3
> Ь1:=1.3;
> z:=0.7;
> t:=time();
ы = 1.3
г = 0.7
г := 3432.022
> E^sum^k/GAMMA^l+a^k^k^^infinity);
> TIME=time()-t;
> ML2(1.3,1.3,0.7);
> TIME=time()-t;
Е- 1.530959770 -0.2083825586 10 *7
TR1E = 464.508
1.707760854
TIME = 0.
Из рассмотренного примера видно, что не только процедура МЬ2 корректно вычисляет функцию типа Миттаг-Леффлера, но и делает это гораздо быстрее, чем по определению (2). Рассмотрим следующий пример.
Пример 3.
> а1:=0.6;
а1 — 0.6
> Ь1:=2;
> z:=0.01;
> t:=time();
Ы - 2
z- 0.01
г:= 3391.102
> E:=sum(z‘k/GAMMA(b1+a1*k),k=0..infinity);
> TIME=time()-t;
> t1:=time();
> ML1(0.6,2,0.01);
1.007036098
> TIME=time()-t1;
В этом примере показано, что численный алгоритм для функции типа Миттаг-Леффлера очень хорошо согласуется с вычислениями по формуле (2).
Заключение
В работе рассмотрен вычислительный алгоритм специальной функции типа Миттаг-Леффлера с параметрами 0 < а < 2 и в > 0. Результаты вычислений имют хорошее согласие с определением функции по формуле (2). Необходимо отметить, что скорость вычисления функции выше, чем скорость вычисления по определению. Это обусловлено тем, что в определении (2) мы грамотно «обрезали» бесконечный ряд погрешностью Р.
Работу можно считать методической, а также использовать как рекомендацию по вычислению функции типа Миттаг-Леффлера в системе MAPLE.
Более детально с алгоритмом вычисления этой функции можно ознакомиться в работе [7], а в работе [8] предложен алгоритм вычисления двухпараметрической функции типа Миттаг-Леффлера.
Библиографический список
1. Нахушев А.М. Дробное исчисление и его применение. М.: Физматлит. 2003. 272 с.
2. Mittag-Leffler G. Sur la repre’sentation analytique d’une branche uniforme d’une fonction monoge‘ne // Acta mathematica. 1905. V. 29. Issue 1. P. 101-181.
3. Джрбашян М.М. Интегральные преобразования и представления функций в комплексной области. М.: Наука, 1966. 671 с.
4. Mainardi F. Fractional relaxation-oscillation and fractional diffusion-wave phenomena chaos, solitons and fractals. 1996. V. 7. № 9. P. 1461-1477.
5. Паровик Р.И. Обобщенное уравнение Матье // Вестник КРАУНЦ. Физ.-мат. науки. 2011. № 2 (3). С. 12-17.
6. Diethelm K., Ford N.J., Freed A.D., Luchko Yu. Algorithms for the fractional calculus: A selection of numerical methods // Comput. Methods Appl. Mech. Eng. 2005. V. 194. P. 743-773.
7. Gorenflo R., Mainardi F., Luchko Yu. Computation of the Mittag-Leffler function £ (z) and its derivatives // Fract. Calculus Appl. Anal. 2002. V. 5. P. 491-518; 2003. V. 6. P. 111-112.
8. Яшагин Н.С. Математическое моделирование и исследование осцилляционных явлений в системах с памятью на основе аппарата дробного интегро-дифференцирования: дис. ... канд. физ.-мат. наук. Самара, 2011. 186 с.
Поступила в редакцию / Original article submitted: 04.10.2012