ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2012, том 52, № 12, с. 2133-2139
УДК 519.644.5
СПЕЦИАЛЬНЫЕ КВАДРАТУРНЫЕ ФОРМУЛЫ ОБРАЩЕНИЯ ПРЕОБРАЗОВАНИЯ ЛАПЛАСА
© 2012 г. А. В. Лебедева, В. М. Рябов
(198504С.-Петербург, Ст. Петергоф, Университетский пр-т, 28, СПбГУ, матем.-механ. ф-т) e-mail: natulich@bk.ru; riabov@VR1871.spb.edu Поступила в редакцию 12.05.2012 г.
Изучаются квадратурные формулы обращения интегрального преобразования Лапласа, приспособленные для обращения изображений, соответствующих длительным медленно изменяющимся оригиналам, характерным для задач линейной вязкоупругости. Доказана сходимость специальных квадратурных формул обращения. Библ. 16.
Ключевые слова: преобразование Лапласа, обращение преобразования Лапласа, квадратурные формулы, сходимость квадратурных формул, задачи линейной вязкоупругости, сингулярные ядра.
1. ВВЕДЕНИЕ
Интегральное преобразование Лапласа Д(р) функции-оригинала/(?):
да
ад = о а, (1)
0
представляет собой мощный инструмент для решения широкого класса прикладных задач математической физики. Одним из его главных достоинств является алгебраизация процедур математического анализа, с помощью которой удается свести интегральные и дифференциальные уравнения к более простым. Кроме того, изображение Лапласа является аналитической функцией в некоторой полуплоскости Яер > у, что позволяет привлечь к исследованию решаемой задачи результаты теории функций комплексного переменного.
Как правило, при решении задач операционными методами наиболее трудным этапом является процесс обращения, т.е. определение оригинала по его изображению. Существуют таблицы соответствия функций-оригиналов и их изображений (см. [1]), теоремы разложения, формула обращения Римана—Меллина, позволяющие теоретически точно находить оригинал. Но решение практических задач часто приводит к изображениям, к которым не могут быть применены эти классические приемы обращения. Следовательно, возникает необходимость разработки и применения приближенных методов.
Наиболее полно возможные подходы к задаче обращения и их реализация описаны в [2]. Обзор других способов обращения, не вошедших в [2], приведен в [3]. Теоретические основы операционного исчисления содержатся в классических работах [1], [4]—[6]. Вопросам приложения операционного исчисления к решению прикладных задач, среди прочих, посвящены фундаментальные труды (см. [7], [8]).
Не существует универсального метода обращения, дающего удовлетворительные результаты для произвольного изображения Д(р). Любой конкретный метод обращения должен учитывать специфику поведения изображения (или функции-оригинала), что прежде всего находит отражение в выборе подходящих систем функций в пространствах оригиналов и изображений, с которыми легко работать и с помощью которых могут быть хорошо приближены заданные образы и оригиналы. Выбор метода обращения существенно зависит от способа задания информации об изображении искомого оригинала. Перечислим типичные ситуации:
1) известны значения изображения Д(р) и его производных в некоторой фиксированной точке, отличной от бесконечности;
2133
2) известны значения изображения F(p) и его производных в некоторой окрестности бесконечно удаленной точки;
3) известны значения изображения Д(р) на вещественной полуоси р > 0;
4) известны значения изображения Д(р) во всей полуплоскости Яер > у.
Выбор подходящих методов обращения для указанных ситуаций, их описание либо отсылка к соответствующей литературе рассмотрены в [3].
Цель настоящей работы состоит в рассмотрении методов обращения с помощью квадратурных формул наивысшей степени точности (КФНСТ) специального вида, обобщающих результаты из [2].
2. КВАДРАТУРНЫЕ ФОРМУЛЫ НАИВЫСШЕЙ СТЕПЕНИ ТОЧНОСТИ
Как известно, обращение преобразования Лапласа задается интегралом Римана—Меллина
С + №
яо = -П. Г )йР, с >у, г > о, (2)
С - №
где у — абсцисса сходимости интеграла Лапласа (1). Напомним, что интеграл (2) понимается в смысле главного значения, он не зависит от с и, в случае разрыва оригинала в точке t, мы получаем полусумму предельных значений оригинала слева и справа от точки t.
Положимр = с + п в формуле обращения (2), тогда ехрр) = ехр(с^ехр(/хО. При фиксированном t первый сомножитель постоянен, а второй пробегает единичную окружность на комплексной плоскости бесконечное число раз. С ростом t первый сомножитель и скорость пробегания окружности вторым сомножителем неограниченно возрастают, так что попытка приблизить интеграл в (2) римановыми суммами вряд ли приведет к цели. Например, в простейшем случае для /(1) = 1 имеем Д(р) = 1/р, так что при любом с > 0 сомножитель ехрС) быстро растет с увеличением t, однако оригинал постоянен и не зависит от с.
Далее для простоты считаем у = 0, чего всегда можно добиться домножением оригинала на соответствующую экспоненту.
Пусть при некотором ж > 0 функция фж(р) = рД(р) регулярна в полуплоскости Яер > 0. Преобразуем интеграл (2) к виду
, c + l да S — 1
f(t) = t— f epp-sф/p/1)dp. (3)
2 n i J
Для приближенного вычисления интеграла (3) выберем попарно различные точкир1, ..., рп в полуплоскости Яер > 0 в качестве узлов квадратурной формулы (КФ) вида
c + l да
—s
П | eppsq(p)dp ~ у Аф(Рк) (4)
2 ni
c — i да k = 1
и потребуем, чтобы формула (4) была точна для функций ф(р) = p-j, j = 0, 2л — 1 , т.е. имела наивысшую степень точности (в таком случае будем говорить, что КФ обладает (2n — 1)-свойством). Это требование равносильно выполнению равенств
л
У Akpk = —Ц, j = 0~2Л—~1. (5)
^ Г(s +j)
k = 1
Такая КФ существует и единственна, все ее узлы pk различны и лежат в правой полуплоскости Rep > 0 (см. [2]).
Заметим, что предельные значения оригинала /(+0), /(+да), если они существуют, могут быть вычислены по формулам
f(+0) = lim pF(p), f(+да) = lim pF(p).
p ^ да p ^ 0
л
По построению формула (4) точна для оригиналов вида I5 1б2и -1(0, где 02п — 1 — любой многочлен степени не выше 2п — 1.
В предположении существования значений /(+0), /(+да) и ограниченности оригинала следует полагать 5 = 1, тогда КФНСТ (4) при I —»- 0 и I —»- да приводит к точному результату, так как сумма коэффициентов КФ (4) равна единице (см. (5)).
К сожалению, КФНСТ (4) плохо приспособлены для обращения изображений, соответствующих медленно протекающим длительным процессам. Так, в задачах линейной вязкоупругости (см. [9]), описывающих напряженное состояние на основе определяющего соотношения Больц-мана—Вольтерра (ниже пространственные координаты для простоты опущены), деформации б и напряжения а связаны соотношением (обобщенный закон Гука)
(
б(t) = 1 a(t) + X \K(t- t)ct(t)dT
E
о
(6)
Первое слагаемое справа в (6) соответствует мгновенной деформации, а второе — наследственной деформации. Как правило, из эксперимента определяется функция ползучести материала — значение правой части (6) при а = const, т.е.
( ' Л
Е( 0 = - ^ 1 + X ¿т). (7)
Важнейшей задачей становится выбор подходящего ядра К интегрального уравнения (6), определяющего функцию ползучести (7). Ядро К должно иметь интегрируемую особенность в точке I = 0. Чаще всего в качестве такового берут дробно-экспоненциальную функцию Работно-ва (см. [9]) (резольвента ядра Абеля)
1 + а к
Эа(в 0 = ?I Г(( 1(!'а)(1+ к)), _1 <а в ° (8)
к = °
Способ определения параметров дробно-экспоненциальной функции по измеренной функции ползучести описан в [10].
Интеграл от этого ядра по полуоси I > 0 должен быть конечным, для чего необходимо р < 0. Не умаляя общности, далее считаем р = —1, и пусть символ Эа(0 означает Эа(—1,1).
В наследственной механике твердого тела наряду с функцией (8) широко используется и интеграл от нее с переменным верхним пределом. Для облегчения использования этих величин составлены таблицы функций (см. [9])
(
а - 1 1%. ^ ч , а + 1
X = I .
F (а, x) = t аЭа(х), F2 (а, x) = t а :|Эа(т) dT, Заметим, что F[(a, x) = E1//a(—x, a), a = a + 1, где
~ п
ЕР(^ = I-+--Г
п=°Г(И + ПР )
есть функция Миттаг—Леффлера (см. [11]).
Однако при решении конкретных задач необходимо вводить в память вычислительной машины части этих таблиц, соответствующие найденным параметрам Эа-функций, которые заранее неизвестны и определяются в процессе решения задачи (и в итоге таковых в таблице может не оказаться). При изменении параметров приходится эту работу проделывать заново, что неудобно и сопряжено с внесением ошибок.
Применяя преобразование Лапласа к уравнению (6), получаем
в(р ) = 1( 1 + т^)^). (9)
Е р -в
В частности, изображение функции ползучести равно
е(р) = Е1 + . (10)
рЕ\ ра +1 - р^
Применяя преобразование Лапласа к уравнению движения среды, получаем второе соотношение между б (р) и а (р), а затем, используя (9), находим б (р) и а (р). Если искомые функции е^) и а(0 ограничены, то можно положить s = 1 и в качестве ф5(р) рассматривать функции ре (р) и ра (р). Они зависят фактически от pa, а = а + 1. Заметим, что для реальных процессов деформирования значение s можно увеличить: так, изображение по Лапласу второго слагаемого в (7), определяющего наследственную деформацию, равно Х/(р(ра — р)), и можно положить s = 1 + а. Искомые решения е(0 и на конечном по t отрезке времени допускают хорошие приближения вида ^ -10(й при 0 < а ^ 1, где — некоторый многочлен, и при уменьшении а скорость их изменения уменьшается.
В таком случае целесообразно вместо КФНСТ построить и использовать обобщенные квадратурные формулы наивысшей степени точности (ОКФНСТ) вида (4), точные для функций
ф(р) = Ра',] = 0, 2п - 1 , или для оригиналов вида ts-102п-1(1а) (02п-1 — произвольный многочлен). Такие формулы были введены в [12] и исследованы в [13]. В частности, в [13] доказана
Теорема 1. Для того чтобы формула (4) была точна для функций ф(р) = р аа,] = 0, 2п - 1 , необходимо и достаточно выполнение двух условий:
1) формула (4) интерполяционная, т.е. точна для функций ф(р) = р ] = 0, п - 1 ;
2) построенный по узлам формулы (4) многочлен
п(*) = П (Х -Рка) (И)
к = 1
удовлетворяет условиям
с +г да
Г Р
| ер
с - г да
е р юп(р )р йр = 0, т = 0, п- 1, с > 0. (12)
В [13] показано, что многочлен (11), удовлетворя
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.