Научная статья на тему 'Особенности вычисления функции типа Миттаг-Леффлера в системе компьютерной математики'

Особенности вычисления функции типа Миттаг-Леффлера в системе компьютерной математики Текст научной статьи по специальности «Математика»

CC BY
643
84
Читать
i Надоели баннеры? Вы всегда можете отключить рекламу.
Ключевые слова
ФРАКТАЛЬНЫЙ ПАРАМЕТРИЧЕСКИЙ ОСЦИЛЛЯТОР / ФУНКЦИЯ ТИПА МИТТАГ-ЛЕФФЛЕРА / КОНТУР ХАНКЕЛЯ / ИНТЕГРАЛЬНОЕ ПРЕДСТАВЛЕНИЕ / FRACTAL PARAMETRIC OSCILLATOR / FUNCTION OF MITTAG-LEFFLER / CIRCUIT HANKEL / INTEGRAL REPRESENTATION

Аннотация научной статьи по математике, автор научной работы — Паровик Роман Иванович

В работе предложены вычислительные алгоритмы для корректного вычисления функции типа Миттаг-Леффлера в математическом пакете MAPLE.

i Надоели баннеры? Вы всегда можете отключить рекламу.

Похожие темы научных работ по математике , автор научной работы — Паровик Роман Иванович

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
Предварительный просмотр
i Надоели баннеры? Вы всегда можете отключить рекламу.

CALCULATION SPECIFIC FUNCTIONS OF MITTAG-LEFFLER IN THE COMPUTER MATHEMATICS

In this work the computational algorithms to correctly calculate the functions of Mittag-Leffler with fractional parameter in mathematical package MAPLE.

Текст научной работы на тему «Особенности вычисления функции типа Миттаг-Леффлера в системе компьютерной математики»

УДК 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] показано, что если выполняются условия

О < а < 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) = < -

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

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

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

> 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

i Надоели баннеры? Вы всегда можете отключить рекламу.