Диссертация (Численное обращение интегрального преобразования Лапласа функций специального вида), страница 8

В этойглаве по порядку рассматриваются возможные пути вычисления, и приводятсяконкретные значения для функций Работнова, вычисленные с помощью каждого метода. Вычисления проводились с помощью математического пакета Maple.Вычисление значений дробно–экспоненциальных функций возможно проводить двумя путями: либо используя явные представления (−1 < α 6 0)Эα (t) = tα∞Xn=0(−1)n tn(1+α),Γ((n + 1)(1 + α))Э∗α (t) =ZtЭα (τ ) dτ,0либо используя их изображенияЭα (p) =1p1+α + 1,∗Эα (p) =1p(pα+1 + 1),применяя методы обращения.Сами функции будем задавать их изображениями по Лапласу, посколькуряд, определяющий дробно–экспоненциальную функцию, плохо сходится ввидурастущих по модулю членов ряда, который к тому же является знакопеременным.

Трудность в использовании изображений функций по Лапласу состоит внахождении оригинала f (t) по изображению F (p), т. е. возникает задача обращения преобразования Лапласа. Есть различные таблицы соответствий междуоригиналами и изображениями, но они охватывают далеко не все встречающиеся на практике случаи или значение оригинала выражается через сложные61функции, которые трудно вычисляемы. В результате точное нахождение оригинала или невозможно, или нецелесообразно.

В связи с этим возникает необходимость в построении приближённых методов обращения преобразования Лапласа, которые были бы применимы для различных случаев.Отдельно можно выделить случай, когда параметр α принимает рациональные значения (α =mn,где 0 <mn< 1 и m и n — взаимно простые числа).В этом случае В.М. Амербаевым были предложены формулы для вычислениядробно–экспоненциальной функции и интеграла от неё [1]. В качестве примераниже приведём формулы для вычисления дробно–экспоненциальной функции.β = −1:a) n — чётное, m — чётноеn−1m(k+1)1 X (−1)k t nЭ mn (−1, t) =+m(k+1)tΓ[]k=0nttZZn−1km(k+1)m(k+1)(−1)1 Xet τ n e−τ dτ + e−t τ n eτ dτ  ++mΓ[ m(k+1) + 1]k=0n−1n0m−222rπ2 X X (−1)k et cos m+m(k+1)m+ 1]r=1 Γ[k=0Ztnτ0m(k+1)ne−τ cos2rπm2rπ4rπdτ,cos+ (t − τ ) sinmm0b) n — чётное, m — нечётноеn−1m(k+1)1 X (−1)k t nЭ mn (−1, t) =+m(k+1)t]Γ[k=0nn−1et X(−1)k+mΓ[ m(k+1) + 1]k=0n−1m−122rπ2 X X (−1)k et cos m+m(k+1)m+ 1]r=1 Γ[k=0nnZtτm(k+1)neZ tτm(k+1)ne−τ dτ +0−τ cos 2rπm 4rπ2rπdτ .cos+ (t − τ ) sinmm0В книге [1] приведены аналогичные формулы для вычисления дробно–экспоненциальной функции интеграла от неё для других случаев (n — нечётное).62Однако стоит заметить, что вычисления с помощью предложенных формул оказываются трудоёмкими, поскольку вычислительные формулы слишкомгромоздки, в них присутствуют интегралы от осциллирующих функций.Теперь рассмотрим способы обращения преобразования Лапласа изображения, зависящего от1paспособами, описанными в Главе 1.633.1.Применение методов деформирования контура к задачеобращения преобразования Лапласадробно–экспоненциальной функции.Рассмотрим применение метода деформирования контура интегрирования описанного в Главе 2 Раздел 2.5 для вычисления дробно–экспоненциальнойфункции.

Изображение Лапласа функции Эα (t) имеет точку ветвления приp = 0, а изображение интеграла от нее еще и простой полюс в той же точке [30]. Для устранения многозначности достаточно выбрать одну из ветвей,что делается стандартным образом: приведенные выше преобразования Лапласа функции Эα (t) и интеграла от нее не имеют особенностей на комплекснойплоскости C \ R− с разрезом вдоль полупрямойR− = {p ∈ C : =(p) = 0, <(p) 6 0}.Пусть в формуле (22) контур интегрирования L состоит из нижнего и верхнего берегов разрезов, соединённых окружностью сколь угодно малого радиусас центром в точке p = 0.Изображения F (p) видаF (p) =1p1+α + 1фактически зависят от pa , так что положим F (p) = Φ1 (pa ) и введем в рассмотрение функцииF ± (t) = Φ1 (ta exp(±iaπ)),t > 0.Очевидно, F + (t) = F − (t) в силу вещественности функции–оригинала.

Воспользуемся полученным в работе [51] следующим результатом:Лемма. Пусть выполнены условия(A) F (p) = o(1) при |p| → ∞, F (p) = o(|p|−1 ) при |p| → 0 равномерно влюбом секторе | arg p| < π − η, π > η > 0;64(Б) существует ε > 0 такое, что для любого ϕ, удовлетворяющего неравенству π − ε < ϕ 6 π, справедливы соотношенияF (r exp(±iϕ))∈ L1 (R+ ),1+r|F (r exp(±iϕ))| 6 α(r),где α(r) не зависит от ϕ и α(r) exp(−δr) ∈ L1 (R+ ) для любого δ > 0.Тогда1f (x) = L (F )(x) =π−1Z∞e−xt = F − (t) dt.(43)0Пусть F (p) = 1/(pa + 1), тогда1ta sin πaF (t) = = a=,t exp(−iπa) + 1 1 + 2ta cos πa + t2a−так что выполнены все условия леммы и формула (43) даетZsin πa ∞ −xtta dtЭα (x) =e=a cos πa + t2aπ1+2tZ0 ∞sinπaz a e−z dza−1=x.πz 2a + 2z a xa cos πa + x2a0(44)При x → 0 последний интеграл в представлении (44) стремится к величинеΓ(1 − a), и с учетом формулы Γ(a)Γ(1 − a) = π/ sin πa при x → 0 из представления (44) получаем Эα (x) ≈ xa−1 /Γ(a), что совпадает с первым членом ряда,определяющего дробно–экспоненциальную функцию.ПоложимxZga (x) =Эα (t) dt.0Изображение этой функции равноGa (p) =1.p(pa + 1)(45)Для нее не выполняется условие (А) леммы (при p → 0 величина |Ga (p)|слишком быстро возрастает).Представим Ga (p) в виде1pa−1Ga (p) = − ap p +1(46)65и положимpa−1Qa (p) = a.p +1(47)Обозначим через qa (x) функцию–оригинал с изображением (47).

Формула(46) означает, что ga (x) = 1 − qa (x).Изображение (47) удовлетворяет условиям леммы, для него находимF − (t) = =ta−1 exp(−iπ(a − 1))ta−1 sin πa=.ta exp(−iπa) + 11 + 2ta cos πa + t2aСледовательно,Zsin πa ∞ −xtta−1 dtqa (x) =e=π1 + 2ta cos πa + t2a0Z ∞sinπaz a−1 e−z dza=x.πz 2a + 2z a xa cos πa + x2a0(48)Из определения (45) следует равенствоga (0) = lim p Ga (p) = 0,p→∞поэтому необходимо qa (0) = 1. Подставив x = 0 в первый интеграл в представлении (48) и сделав замену ta = z, придем к табличному легко вычисляемомуинтегралу и таким образом убедимся в справедливости равенства qa (0) = 1 привсех a > 0.Итак, наши задачи обращения преобразования Лапласа свелись к вычислению интеграловZsin πaЭα (x) = xa−1πsin πaqa (x) = xaπ∞0Z0∞z a e−z dz,z 2a + 2z a xa cos πa + x2az a−1 e−z dz.z 2a + 2z a xa cos πa + x2a(49)(50)Для их приближенного вычисления можно применить квадратурные формулы типа Гаусса [25] с весом Лагерра z a e−z для первого интеграла и с весомz a−1 e−z для второго интеграла.

Однако с уменьшением a точность формул будетуменьшаться.66Поэтому для приближенного вычисления интегралов (49) и (50) построимобобщенные квадратурные формулы видаZ ∞nXβ −zAk f (zk ),z e f (z) dz ≈0β > −1,(51)k=1точные для функций f (z) = z am , m = 0, 1 . . . , 2n − 1.В работе [27] доказана следующаяТеорема 9. Для того чтобы формула (51) была точна для функцийf (z) = z am , m = 0, 1 . . . , 2n − 1, необходимо и достаточно выполнение двухусловий:1) формула (51) интерполяционная;2) построенный по узлам формулы (51) многочленωn (z) =nY(z − zka )(52)k=1удовлетворяет условиямZ ∞z β e−z ωn (z a )z am dz = 0,m = 0, 1, .

. . , n − 1.(53)0Д о к а з а т е л ь с т в о этой теоремы проводится точно так же, как в случае классических формул типа Гаусса [25].Покажем, что многочлен (52), удовлетворяющий условиям (53), существует и определяется однозначно.После замены переменной z a = x условия (53) принимают видZ ∞x(β+1)/a−1 exp(−x1/a ) ωn (x) xm dx = 0, m = 0, 1, . . . , n − 1.0Функция w(x) = x(β+1)/a−1 exp(−x1/a ) обладает свойствами веса на полуоси(0, ∞), поскольку (β + 1)/a > 0, следовательно, искомый многочлен существует и единствен, а его корни, т.

е. zka , k = 1, 2, . . . , n, попарно различны и положительны. Все коэффициенты формулы положительны. Итак, квадратурнаяформула типа Гаусса вида (51) существует.67Опишем способ вычисления узлов и коэффициентов формулы (51). Будемискать многочлен (52) в видеωn (z) = z n + b1 z n−1 + · · · + bnс неизвестными коэффициентами bk .Условия (53) приводят к системе линейных алгебраических уравненийnXΓ(β + (k + n − j)a + 1) bj = −Γ(β + (k + n)a + 1),k = 0, 1, . .

. , n − 1.j=1Ее решение существует и единственно, как показано выше. Далее находим корни уравнения ωn (z) = 0, т. е. числа zka , k = 1, 2, . . . , n. Коэффициенты формулы(51) определяем из системы уравненийnXAk (zka )j−1 = Γ(β + (j − 1) a + 1),j = 1, 2, . . . , n.k=1Заметим, что узлы и коэффициенты формулы (51) вещественны.Описанный метод был реализован в виде программы в математическомпакете Maple. Текст программы приведён в Приложении Программа 7. Здесьприведем результаты вычисления функций F1 (α, x), F2 (α, x).Результаты вычисления F1 (α, x) при α = −0.4.Таблица 17nxF1F1табл.300.15.724888 · 10−15.725716 · 10−130 0.155.298362 · 10−15.298110 · 10−1300.24.9090271 · 10−1 4.909080 · 10−1300.82.167449 · 10−13012.167449 · 10−11.7110228 · 10−1 1.711022 · 10−168Результаты вычисления F1 (α, x) при α = −0.65.Таблица 18nxF1F1табл.3019.124 · 10−29.124 · 10−22041.2219 · 10−2 1.222 · 10−23041.2219 · 10−2 1.222 · 10−230 400 1.5759 · 10−6 1.576 · 10−6Результаты вычисления F2 (α, x) при α = −0.4.Таблица 19nxF2F2табл.30 0.2 9.5903 · 10−1 9.5907 · 10−130 0.8 6.5302 · 10−1 6.5302 · 10−1305.8667 · 10−1 5.8667 · 10−11Результаты вычисления F2 (α, x) при α = −0.65.Таблица 20nxF2F2табл.3015.507 · 10−1 5.507 · 10−12042.101 · 10−1 2.101 · 10−13042.101 · 10−1 2.101 · 10−120 400 2.495 · 10−3 2.495 · 10−3Как видно из результатов, приведённых в таблицах 17–20, наиболее точные значения (по сравнению с табличными) достигаются при значениях аргумента x > 1.693.2.Обращение преобразования Лапласа по значениямизображения на вещественной осиВ работах [29] и [42] предложен один из вариантов обращения преобразования Лапласа по значениям изображения на вещественной оси, в основе котороголежит следующая теорема:Теорема 10 ([63], с.

385). Пусть f (t) ∈ L(0, ∞) и ее преобразованиеЛапласа равно∞ZF (p) =exp(−p t)f (t) dt.0Тогда для почти всех положительных значений t имеет место формула обращенияZ(−1)n tn−1 ∞ −p t 2n−1 (n)f (t) = lime pF (p) dp.(54)n→∞ n!(n − 2)! 0Возможность применения этой формулы ограничена в связи с тем, что вней присутствует производная изображения. Однако интегрированием по частям она может быть представлена иначе [29]:Z ∞f (t) = lime−p t P2n−1 (p t)F (p) dp,n→∞где0n (−1)n−1 (2n − 1)! X n (−p)2n−j−1.P2n−1 (p) =n!(n − 2)!j(2n−j−1)!j=0Как показано в работе [29], этот многочлен выражается через многочлены Лагерра общего видаLn (x ; α) =nXmn + α xmn − m m!(−1)m=0(55)простой формулой:P2n−1 (p) =pn−1Ln (p ; n − 1).(n − 2)!В результате этих преобразований утверждение (54) теоремы запишетсяв видеf (t) = lim fn (t),n→∞1fn (t) =t(n − 2)!Z0∞e−p pn−1 Ln (p ; n − 1)F (p/t) dp.70Для приближенного вычисления последнего интеграла применим квадратурную формулу типа Гаусса с весом e−p видаZ∞−pe ψ(p) dp ≈0mXAk ψ(pk ),(56)k=1точную для всех многочленов степени не выше 2m − 1.Такой метод обращения преобразования Лапласа, пригодный и для определения величин скачков оригинала в точках разрыва первого рода, был предложен и исследован в работе [37].