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

. . , k. Положим mj =11++,n1 nkдостигает1xj ,максимума,т. е.j = 1, . . . , k, а затем поmjm1 ,j = 1, . . . , k, построим числа cjk вида (18). Тогда, как покаkPзано в работе [32], величина|cjk | принимает минимальное значение, равноеj=1(−1)k−1 Tek−1 (0).Числа mj не целые, и поэтому следует в качестве искомых номеров приближений по Виддеру взять ближайшие целые к ним, т.

е. положить nj =bmj + 12 c,j = 1, 2, . . . , k. Такие номера будем называть чебышевскими.Необходимость вычислять оптимальные значения для параметров r, m какэто описано выше, увеличивает трудоёмкость применения данного метода. Оптимальным значением параметров r, m являются r = 0.95, m = 950, параметрk — число натуральных чисел n1 , n2 , . . . , nk , принимался равным 3. При этом с33увеличением m, k, а также длины промежутка nk − n1 точность вычислений неувеличивалась.Метод, описанный в этой главе, был применён для вычисления функцииползучести (пространственная переменная x опущена)Ztε(t) = 1 + 4Э−0.8 (τ )dτ0по её изображениюε(p) =41+.p p(p0.2 + 1)Программа для вычисления представлена в Приложении Программа 3. Результаты вычислений находятся в Таблицах 8–10.Результаты вычисления ε(t) при α = −0.8,t = 0.001.Таблица 8nε(t)ε(t)точное100 1.866991.86696150 1.866961.86696200 1.866961.86696Результаты вычисления ε(t) при α = −0.8,Таблица 9nε(t)ε(t)точное100 3.115653.11559150 3.115603.11559200 3.115593.11559t = 1.34Результаты вычисления ε(t) при α = −0.8,t = 106 .Таблица 10nε(t)ε(t)точное100 4.793554.79352150 4.793474.79352200 4.793474.79352Применяя метод Виддера к задаче нахождения функции–оригинала поизображению по Лапласу, можно добиться более точных значений по сравнениюс рассмотренными в предыдущих параграфах методами (КФНСТ, ОКФНСТ).Но, чтобы добиться точности вычислений, необходимо увеличивать точностьпромежуточных вычислений — это проявление неустойчивости исходной задачи, и если её наращивать в КФ, то они также будут давать хорошие результаты.Метод Виддера для обращения интегрального преобразования Лапласарассматривался и изучался в работах [43], [27].351.5.Деформирование контура интегрированияРассмотрим обращение преобразования Лапласа, исходя из формулы Римана–Меллина1f (t) =2πiZc+i∞ez t F (z) dz,c > γ,(21)c−i∞где γ — абсцисса сходимости интеграла Лапласа.

Напомним, что интеграл (21)понимается в смысле главного значения, он не зависит от c и в случае разрываоригинала в точке t мы получаем полусумму предельных значений оригиналаслева и справа от точки t.Положим в формуле обращения (21) z = c + iτ, тогда ez t = ect eiτ t . Прификсированном t первый сомножитель постоянен, а второй пробегает единичную окружность на комплексной плоскости бесконечное число раз. С ростомt первый сомножитель и скорость пробегания окружности вторым сомножителем неограниченно возрастают, поэтому римановыми суммами приблизитьинтеграл в (21) скорее всего не получится.Для преодоления этих трудностей заменим линию интегрирования в (21)эквивалентным контуром L = {z | z = l(u), u ∈ (−∞, +∞)}, который начинается и заканчивается в левой полуплоскости так, что <(z) → −∞ на обоих егоконцах. Такая замена возможна при выполнении условий:1) внутри контура L содержатся все особенности изображения F (z);2) |F (z)| → 0 равномерно в полуплоскости <(z) 6 γ при |z| → ∞.

Далеевсюду предполагается, что эти условия выполняются.Запишем интеграл (21) в виде1f (t) =2πiZezt F (z) dz,t > 0,(22)Lа затем положим z = l(u), в результате формула (22) принимает видZ ∞f (t) =Gt (u) du,−∞(23)36где1 tl(u) 0e l (u)F (l(u)).(24)2πiФункция (24) не имеет особенностей как на линии интегрирования, такGt (u) =и в некоторой “полосе”, содержащей внутри себя линию интегрирования.

Длявычисления интеграла (23) воспользуемся формулой трапеций с бесконечнымчислом узлом (см. [45], [49]). Заметим, что в (23) интегрирование происходит побесконечному промежутку, а формула погрешности аппроксимации для классической формулы трапеции применима лишь для интегрирования по конечномупромежутку.Так, классическая составная формула трапеций выглядит следующим образом (составная формула трапеций Котеса):!Z bn−1Xf0 + fnf (x) dx = hfi + En (f ),+2ai=1где h — шаг сетки, En (f ) — погрешность аппроксимации задаётся выражением00f (ξ)En (f ) = −(b − a)h2 ,12где ξ ∈ (a,b).Однако, в работе [45] показано, что формула трапеций может давать хорошуюточность и для нашего случая (промежуток интегрирования — бесконечный),при этом скорость сходимости зависит от ширины области регулярности функции Gt (u) и шага численного интегрирования.Задача вычисления оценок погрешности в результате применения формулы трапеций с бесконечным числом узлов рассмотрена в работе [62].Пусть w = u + iv, u, v ∈ R, функция g(w) аналитична в полосе −d 6 v 6 cдля некоторых c > 0, d > 0, и равномерно в этой полосе g(w) → 0 при |w| → ∞c такой скоростью, что существует интегралZ +∞|g(u + ir)|du, r ∈ [−d, c].−∞Для вычисления интегралаZI=+∞g(u)du−∞37применим формулу трапеций c бесконечным числом узлов:I ≈ Ih = h+∞Xg(kh),k=−∞и с конечным числом (2N + 1) узлов:I ≈ Ih ,N = hNXg(kh).k=−NПоложим DE = |I − Ih |, T E = |Ih − Ih ,N |.Теорема 8.

Пусть для некоторых положительных чисел M+ , M− справедливы неравенства+∞Z|g(u + ir)|du 6 M− ,r ∈ [−d, 0],−∞Z+∞|g(u + is)|du 6 M+ ,s ∈ [0, c].−∞Тогда|I − Ih | 6 DE+ + DE− ,гдеDE+ =M+,e2cπ/h − 1DE− =M−.e2dπ/h − 1Доказательство теоремы 8 приведено в книге [42].Константы M+ , M− не зависят от шага h численного интегрирования изначения N , и потому в дальнейшем ограничимся качественным поведениемоценок с помощью соотношенийDE+ = O(e−2cπ/h ), DE− = O(e−2dπ/h ),Для погрешности T E при постоянном h полагаемT E = O(|g(hN )|),N → ∞.h → 0.38Далее для избранного контура интегрирования приравниваем характеристики величин DE+ , DE− , T E, что приведет к некоторому способу выбора параметров контура и шага интегрирования в зависимости от N и возможностиполучения полной погрешности метода.Результаты, описанные выше, применим для вычисления интеграла (23),положив g(w) = Gt (w).В работах [40] и [61] были предложены две линии интегрирования дляинтеграла (22), однако для них полоса регулярности имеет переменную ширину и стремится к нулю при неограниченном возрастании модуля переменнойинтегрирования.В работе [62] предложены другие, более эффективные контуры, для которых получены априорные оценки погрешности.

В работе [28] исследовалисьпараболический и гиперболический контуры из [62] и даны оценки погрешностидля каждого случая.Так, в работе [62] в качестве контура предлагается брать кривуюz(u) = µ(1 − u2 ) + 2iµu,(25)где µ > 0, тогда формула обращения примет вид:Z ∞1f (t) =ez(u)t F (z(u))z 0 (u) du.2πi −∞Прямые w = u − id при v = d,d > 0, и w = u + ic при v = c,0<c<1при отображении (25) переходят в параболы:z(u) = µ((1 + d)2 − u2 ) + 2iuµ(1 + d),z(u) = µ((1 − c)2 − u2 ) + 2iuµ(1 − c).Предположим, что все особые точки изображения расположены в конечной части полуплоскости <(z) < 0, и не все они вещественны.Существует оптимальное значение параметра d, при котором достигаетсянаилучшая оценка [62] 2− thπ2 µ + 2πhDE− = O e,h → 0.39В нашем предположении ветви параболы не смыкаются, таким образом,наибольшее допустимое значение c удовлетворяет неравенству 0 < c < 1.Приравняв показатели величин DE− , DE+ , T E, найдем h, µ:p1 + 4c(c + 1)ππNph=, µ==.N2(c + 1)th 2(c + 1) 1 + 4c(c + 1) tНапомним, что число c здесь не произвольно.

Например, в случае f (t) =sin(ωt) изображение F (z) = ω/(z 2 + ω 2 ) имеет особые точки ±iω, которыедолжны лежать между точками пересечения параболы z(u) = µ((1 − c)2 − u2 ) +2iuµ(1−c) с мнимой осью, т. е. должно выполняться неравенство 2µ(1−c)2 > |ω|или равносильное ему неравенствоp|ω| t(c + 1) 1 + 4c(c + 1)N>.π(1 − c)2Погрешность метода в таком случае есть величина!!2πc2π= O exp − pN,O exp − ch1 + 4c(c + 1)N → ∞.(26)В качестве гиперболического контура предлагается брать кривуюz(u) = µ(1 + sin(iw − α)), w = u + iv.Аналогично случаю параболического контура имеем оптимальное значение параметра d, при котором достигается наилучшая оценка [62]DE− = O (exp(µt − 2πα/h)) ,h → 0,T E = O (exp(µt(1 − sin(α) cosh(hN )))) ,N → ∞.Приравняв показатели величин DE− , DE+ , T E, найдем h, µ:A(α),h=N4πα − π 2 Nµ=,A(α) tгдеA(α) = cosh−12α.(4α − π) sin(α)40Заметим, что в описанных выше случаях получен лишь порядок оценкиубывания погрешности.В статье [21] приводится кусочно прямолинейный контур (см.

Рис. 1,Рис. 2), при использовании которого получены явные оценки погрешности. Далее для более полного понимания сути работы [21] и ввиду важности полученных в статье результатов приведём здесь её основные результаты.Рис. 1. Первый контурРис. 2. Второй контурТак, в [21] предлагается рассмотреть кусочно прямолинейный контур, состоящий из трёх участков.

Обозначим их L1 , L2 , L3 . На участке L2 параметр zдля обоих контуров изменяется по закону z = xbi, x ∈ [−1, 1], b > 0. На участкеL1 имеем z = x − bi, x ∈ (−∞, 0], а на участке L3 z = −x + bi, x ∈ [0, ∞).Заметим, что первый контур является частным случаем второго (еслиуглы наклона участков L1 , L3 второго контура равны нулю). Выбор контураиз предложенных в [21] зависит от расположения особых точек изображения,однако характер исследования одинаков, и потому для простоты ограничимсярассмотрением первого контура.Предполагается, что все особые точки изображения расположены в левой полуплоскости, поэтому в формуле обращения Римана–Меллина (21) можно положить c = 0.

В противном случае вместо F (z) рассмотрим функциюF (z + a), a > 0 (она соответствует функции–оригиналу f (t) exp(−at)) с особыми точками в левой полуплоскости.41Используя предложенный контур для вычисления (22), получимZZZZ 11f (t) =ez t F (z) dz =++ez t F (z) dz.2πi L2πiL1L2L3После замены переменной z получимZZ ∞zt−ibte F (z) dz =ee−xt F (−x − bi) dx,L0Z 1Z 1ez t F (z) dz =bieibtx F (ibx) dx,−1ZL2Z ∞ztibte F (z) dz = − ee−xt F (−x + bi) dx.L30В работе [21] дано подробное описание способа вычисления интегралов,входящих в приведённые выше формулы.