Применение вариантов шагового метода численного обращения преобразования Лапласа Текст научной статьи по специальности «Математика»

Механика

Вестник Нижегородского университета им. Н.И. Лобачевского, 2014, № 1 (2), с. 249-255

УДК 539.3+519.6

ПРИМЕНЕНИЕ ВАРИАНТОВ ШАГОВОГО МЕТОДА ЧИСЛЕННОГО ОБРАЩЕНИЯ ПРЕОБРАЗОВАНИЯ ЛАПЛАСА

© 2014 г. Л.А. Игумнов, Я.Ю. Ратаушко

НИИ механики ННГУ им. Н.И. Лобачевского

igumnov@mech.unn.ru

Поступила в редакцию 25.09.2013

Рассматривается задача построения шаговых по времени методов численного обращения преобразования Лапласа, основанных на теореме операционного исчисления об интегрировании оригинала. Шаговая схема определяется выбором квадратурной формулы и схемой численного решения задачи Коши, порождаемой интегралом, возникающим в рамках рассматриваемого подхода.

Ключевые слова: обращение преобразования Лапласа, шаговый метод, переменный шаг, схема Рун-ге-Кутты.

Введение

Оригинальный подход к построению шаговых по времени схем метода граничных элементов описан в работах [1-3]. Ключевой проблемой построения гранично-элементной схемы [3-5] является проблема численного обращения интегрального преобразования Лапласа. В работах [6-8] применён шаговый метод численного обращения преобразования Лапласа. В [9] приводится модификация метода, позволяющая учитывать специфику интегрирования сильно осциллирующих функций. Расширение подхода дают также работы [10, 11], использующие схемы Рунге-Кутты для решения задачи Коши, порождённой специальной процедурой обращения преобразования Лапласа. Работа посвящена вопросу рассмотрения и сравнения модификаций шагового метода, основанных на изложенных в [9-11] идеях.

Постановка задачи и метод решения

Прямое и обратное интегральные преобразования Лапласа соответственно определяются формулами:

ад

У(я) = 1 у(г)е"'—г,

0

у(г) =- | у^У^я,

2п ю-/ад

где я = ю + /ю - комплексная переменная, введённая в полуплоскости ю > ю0 .

Рассмотрим метод, опирающийся на теорему об интегрировании оригинала - шаговый метод

численного обращения преобразования Лапласа.

Пусть

у(г) = } / (т)—т.

(1)

Заменим интеграл (1) квадратурной суммой, весовые множители которой определяются с помощью изображения по Лапласу / и линейного многошагового метода (изложение идёт с учётом результатов [1-3]):

п

у(0) = 0, у(пМ) = (ДО, п = 1.....N, (2)

к=1

где

п-п Ь-1 _ Ю п (Д) = — I /

ь 1=0

(

2л Л

/1

у(Яе ь ) Д

(3)

Аппроксимация, используемая при выводе формул (2), (3), основана на применении линейного многошагового метода (с характеристической функцией у(г)) для решения возникающей в процессе преобразования интеграла (1) задачи Коши для дифференциального уравнения первого порядка:

— х(г) = ях (г) + С, х(0) = 0.

—

(4)

На выбор многошагового метода накладываются следующие условия: метод должен быть порядка точности р > 1, являясь строго нуль-

устойчивым или Л-устойчивым. Функция /(я) должна быть ограничена в правой полуплоскости относительно прямой (с - /ад, с + /ад), то есть:

\/(я) < К|я|при К < ад, ц > 0.

0

ь

е

Если функция f (s) аналитична и ограниче-

на в области

|arg(s - с) < п - ф,

где ф < ^, критерий устойчивости может быть

ослаблен до A(a) -устойчивости.

К соответствующим примерам многошаговых методов относятся методы дифференцирования назад порядка p < 6: для Л-устойчивого метода дифференцирования назад второго порядка (а = 90°) можем записать: у(г) = 3/2 - 2z + z2 /2.

При условии того, что функция f в уравнении (3) вычисляется с некоторой погрешностью е, выбор L = N и Rn = 4е допускает погрешность вычисления ш п порядка

Модификация шагового метода с переменным шагом

Модификация формулы (3) для вычисления шп с переменным шагом и линейной аппрокси-

Модификация шагового метода на узлах схемы Рунге-Кутты

Ещё один возможный вариант модификации традиционного метода с целью повышения качества результата при сохранении затрат - сгущение расчётных узлов в пространстве частот путём замены многошагового метода, использованного для решения задачи Коши (4).

Рассмотрим метод Рунге-Кутты, записанный с помощью таблицы Бутчера:

е\АТ

-1—, Л е Rmxm , Ь,С еЯп. ЬТ

Для корректной формулировки шаговой схемы должны быть выполнены следующие условия [11]:

1. Метод Рунге-Кутты должен быть А-устойчивым;

2. |Л(г)| < 1 при 7 ф 0 , где R(z) = 1 +

Т 1 т

+ zb1 (I - гЛ)-1(1,...,1)т - функция устойчиво-

сти;

3. R(<x>) = 0;

4. 3AT1.

Для простоты

будем рассматривать

мацией функции выглядит следующим образом: b A = (0, -,0,1), тогда метод автоматически L-

Ю (At) =

R

-n L-1

2 X f

2п к=0

fy(Re )>

+ f

y( Re

m<Pk+i)

At

-m^k+1

e +

At

/

Фк+1 - Фк (5)

2

Для случаев, когда }(у(Reгnф)/ А/)/ еп -сильно осциллирующая функция, в сочетании с (5) целесообразно использовать комбинированную формулу, учитывающую специфику интегрирования таких функций [9]:

устойчив.

Внедряя формализованный метод Рунге-Кутты вместо линейного многошагового метода для решения задачи Коши дифференциального уравнения, получим [11]:

п

Уо = 0, Уп = Ътл~х (А/), п = 1,...,Ы, к=1

R-n L-1 -®n (At) = — X f

L 1=0

(

:,2Л Л

A(Re L ) At

-ыЦ-

•¡(At)=^ X фк+1 -фк

2п к =0

2

A(w)f

у( Re'nn ) At

+ D2(w)f

..Фк+Фк+1

rj(Rein,fk+1)^

At

(6)

Фк+1 -Фк 2

sin w wcosw - sin w | | -+---i при I w > w2,

D 2 (w) = <j w w

+wi I I ^

e при I w < w2.

Формулы (5), (6) позволяют перераспределять расчётные узлы по промежутку изменения ф для получения большей точности результата.

А(г) = Л'1 - гЛ^ЩЪ Л-.

В качестве схем метода Рунге-Кутты, удовлетворяющих сформулированным условиям, выберем схемы Радо и Лобатто [11].

Численные результаты

Рассматривается решение задачи о действии осевой силы Е=1Н/м на пороупругий стержень [3]. Отклики перемещений и давлений, вызванные силой Е, наблюдаются в точке у = 3 м стержня длиной I = 3 м . Численные результаты получаются для материала с параметрами: К = 4.8-109 Н/м2 , О = 7.2 • 109 Н/м2,

e

e

2

Ю

X

X

ря = 2458кг/м3, ф = 0.19,К„ = 3.6-1010 Н/м2 , рг = 1000кг/м3, Кг = 3.3-109 Н/м2, к =

= 1.9 -10-10 м4/(Н • с). Здесь К, О - константы упругости, р я, р у - плотности упругого скелета и наполнителя, ф - пористость, к - проницаемость, , К у - объемные модули упругого

скелета и наполнителя. Выберем отрезок времени 0.02 с, что составляет приблизительно 5.5 периодов функций по времени. При расчётах использован коэффициент Я = 0.997.

На рис. 1 , 2 приведён вид действительной и мнимой частей спектра перемещений и давлений соответственно.

Обозначим через Ь общее число узлов (включая дополнительные для схем Рунге-Кутты) по углу, через N - по времени. Рассмот-

рим равномерное разбиение промежутка изменения ф для модификации шагового метода, основанной на схемах Лобатто и Радо, и кусочно-равномерную сетку на промежутках [0, л /2], [л / 2,3л / 2], [3л / 2,2л] для модификации с переменным шагом, основанной на методе Эйлера (2), (5). Результаты обращения для отклика перемещений при N=¿/2=600 представлены на рис. 3. Красная кривая соответствует схеме с переменным шагом, количество узлов на промежутках - 560, 80, 560; синяя - схеме Лобатто; зелёная - схеме Радо. Здесь и далее для сравнения чёрным цветом приведён более точный вид решения. Аналогично на рис. 4 показаны результаты для N=¿/2=1200, на рис. 5 -для N=¿/2=2400, с сохранением пропорции количества узлов на промежутках для схемы с переменным шагом. При измельчении сеток на

-2.x Ю"10

Рис. 3

-5.x М-11

-1.x Ю10

-их Ю"10

-2.5 х 10-">

11-1

0.020

Рис. 4

-2.x Ю-10

-2.5 х Ю"10

'Л

-1.Х10-"-

-2.x Ю'10-

Рис. 5

и

и

и

промежутках до 2240, 320, 2240 узлов на графике появляются характерные осцилляции, поэтому для подсчётов использована модификация на основе формул интегрирования сильно осциллирующих функций (2), (6).

На вставке в рис. 5 в конце расчётного промежутка видна потеря устойчивости схем, основанных на методах семейства Рунге-Кутты, поэтому на рис. 6а,б дополнительно приведён вид двух предпоследних пиков решения.

0.0142 0.01-- 0.0146 0.0143 0.0150

0.0136 0.0155 0.01« 0.0162 0.0164

-8.x 10"

-1.x

-1.5 х 10"™-

-2.4 х 10"1с-

-2.6х 10-1:'-

а)

б)

Рис. 6

-0.4 -

-0.4 -

Рис. 8

u

P

P

На рис. 7, 8 представлены полученные графики давления при N=¿/2=600, N=¿/2=1200 соответственно; цветовое обозначение использованных схем сохранено.

Схема Радо независимо от числа узлов аппроксимации даёт наименьший сдвиг по времени (запаздывание или опережение) и в целом в рассмотренных условиях незначительно уступает схеме с переменным шагом, основанной на методе Эйлера; схема Лобатто даёт худшие результаты. На малом числе узлов подход на основе методов семейства Рунге-Кутты отстаёт в точности аппроксимации от подхода с переменным шагом, однако схема Радо лучшим образом передаёт качественный характер решения. С дальнейшим измельчением сеток схемы на основе методов Рунге-Кутты предоставляют немного лучшие результаты, но первыми теряют устойчивость при кратном увеличении N и ¿.

Список литературы

1. Lubich C. Convolution Quadrature and Discre-tized Operational Calculus. I // Numerische Mathematik. 1988. № 52. P. 129-145.

2. Lubich C. Convolution Quadrature and Discre-tized Operational Calculus. II // Numerische Mathematik. 1988. № 52. P. 413-425.

3. Schanz M. Wave Propagation in Viscoelastic and Poroelastic Continua. Berlin Springer, 2001. 170 p.

4. Баженов В.Г., Игумнов Л.А. Методы граничных интегральных уравнений и граничных элементов в решении задач трехмерной динамической тео-

рии упругости с сопряженными полями. М.: Физмат-лит, 2008. 352 с.

5. Угодчиков А.Г., Хуторянский Н.М. Метод граничных элементов в механике деформируемого твердого тела. Казань: Изд-во Казанского ун-та, 1986. 295 с.

6. Белов А.А., Игумнов Л.А., Литвинчук С.Ю. Развитие метода граничных элементов для решения трехмерных контактных нестационарных динамических задач теории упругости // Проблемы прочности и пластичности: Межвуз. сборник. Н.Новгород: Изд-во ННГУ, 2007. Вып. 69. С. 125-136.

7. Аменицкий А.В., Белов А.А., Игумнов Л.А., Литвинчук С.Ю. Гранично-элементное моделирование на основе квадратур сверток динамического состояния составных упругих тел // Вычислительная механика сплошных сред. Пермь: Изд-во ИМСС УрО РАН. 2008. Т. 1. № 3. С. 5-14.

8. Аменицкий А.В. Развитие методов граничных элементов для численного моделирования динамики трехмерных однородных пороупругих тел: Автореф. дис... канд. физ.-мат. наук. Н.Новгород, 2010. 20 с.

9. Белов А.А., Игумнов Л.А., Литвинчук С.Ю. Гранично-элементная методика на основе модифицированного метода квадратур сверток в динамических задачах упругих тел // Проблемы прочности и пластичности: Межвуз. сборник. Н.Новгород: Изд-во ННГУ, 2008. Вып. 70. C. 150-158.

10. Banjai L. Multistep and multistage convolution quadrature for the wave equation: Algorithms and experiments // SIAM J. Sci. Comput. 2010. № 32. P. 29642994.

11. Banjai L., Messner M., Schanz M. Runge-Kutta convolution quadrature for the boundary element method // Comput. Methods Appl. Mech. Engrg. 2012. P. 90-101.

APPLICATION OF TIME-STEP METHOD VARIANTS FOR NUMERICAL INVERSION

OF THE LAPLACE TRANSFORM

L.A. Igumnov, Ya. Yu. Rataushko

The problem of constructing time-step methods for numerical inversion of the Laplace transform, based on the original function integration theorem, is considered. The stepping scheme is determined by the choice of the quadrature formula and the numerical solution scheme for the Cauchy problem, which arises for the integral appearing in the framework of the observed approach.

Keywords: Laplace transform inversion, time-step method, variable step, Runge-Kutta scheme.

References

1. Lubich C. Convolution Quadrature and Discre-tized Operational Calculus. I // Numerische Mathematik. 1988. № 52. P. 129-145.

2. Lubich C. Convolution Quadrature and Discre-tized Operational Calculus. II // Numerische Mathematik. 1988. № 52. P. 413-425.

3. Schanz M. Wave Propagation in Viscoelastic and Poroelastic Continua. Berlin Springer, 2001. 170 p.

4. Bazhenov V.G., Igumnov L.A. Metody granich-nyh integral'nyh uravnenij i granichnyh jelementov v

reshenii zadach trehmernoj dinamicheskoj teorii upru-gosti s soprjazhennymi poljami. M.: Fizmatlit, 2008. 352 s.

5. Ugodchikov A.G., Hutorjanskij N.M. Metod gra-nichnyh jelementov v mehanike deformiruemogo tver-dogo tela. Kazan': Izd-vo Kazanskogo un-ta, 1986. 295 s.

6. Belov A.A., Igumnov L.A., Litvinchuk S.Ju. Razvitie metoda granichnyh jelementov dlja reshenija trehmernyh kontaktnyh nestacionarnyh dinamicheskih zadach teorii uprugosti // Problemy prochnosti i plastich-nosti: Mezhvuz. sbornik. N.Novgorod: Izd-vo NNGU, 2007. Vyp. 69. S. 125-136.

7. Amenickij A.V., Belov A.A., Igumnov L.A., Litvinchuk S.Ju. Granichno-jelementnoe modelirovanie na osnove kvadratur svertok dinamicheskogo sostojanija sostavnyh uprugih tel // Vychislitel'naja mehanika sploshnyh sred. Perm': Izd-vo IMSS UrO RAN. 2008. T. 1. № 3. S. 5-14.

8. Amenickij A.V. Razvitie metodov granichnyh jelementov dlja chislennogo modelirovanija dinamiki trehmernyh odnorodnyh porouprugih tel: Avtoref. dis... kand. fiz.-mat. nauk. N.Novgorod, 2010. 20 s.

9. Belov A.A., Igumnov L.A., Litvinchuk S.Ju. Granichno-jelementnaja metodika na osnove modificiro-

vannogo metoda kvadratur svertok v dinamicheskih zadachah uprugih tel // Problemy prochnosti i plastich-nosti: Mezhvuz. sbornik. N.Novgorod: Izd-vo NNGU, 2008. Vyp. 70. C. 150-158.

10. Banjai L. Multistep and multistage convolution quadrature for the wave equation: Algorithms and experiments // SIAM J. Sci. Comput. 2010. № 32. P. 29642994.

11. Banjai L., Messner M., Schanz M. Runge-Kutta convolution quadrature for the boundary element method // Comput. Methods Appl. Mech. Engrg. 2012. P. 90101.