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

А также приведены оценки погрешности для вычисления каждого из интегралов.В случае интегралов видаZ1eiωx g(x) dx(27)−1предлагается производить вычисления следующим образом.Воспользуемся предложенным в работе [3] способом приближенного вычисления интегралов вида (27): дана квадратурная формула Гаусса [15] с nузлами алгебраической степени точности 2n − 1Z 1nXv(x) dx ≈Aj v(xj ).1(28)j=1Её узлы суть корни многочлена Лежандра Pn (x).Построим интерполяционный многочлен Ln−1 (x) для функции g(x) позначениям g(x1 ), . .

. , g(xn ) и представим его в виде разложения по многочленамЛежандраLn−1 (x) =n−1Xck Pk (x),k=0коэффициенты которого точно вычисляются с помощью формулы Гаусса (28):n2k + 1 XAj g(xj )Pk (xj ).ck =2 j=142Погрешность интерполяции представима в видеg (n) (ξ) 2n (n!)2g(x) − Ln−1 (x) =Pn (x), ξ ∈ [−1, 1].n!(2n)!(29)Полагая в формуле (27) g(x) ≈ Ln−1 (x) и используя известное равенствоrZ 12πJk+1/2 (ω),eiωx Pk (x) dx = ikω−1где Jν (ω) — функции Бесселя, придем к интерполяционной квадратурной формулеZ1e−1iωxg(x) dx ≈nXrDj (ω)g(xj ), Dj (ω) = Ajj=1n−12π X 2k + 1 ki Jk+1/2 (ω)Pk (xj ),ω2k=0(30)погрешность которой равнаZ1Rn (g) =eiωx (g(x) − Ln−1 (x)) dx−1и в силу преставления (29) оценивается неравенством|Rn (g)| 621√max |g (n) (x)|.(2n − 1)!! 2n + 1 [−1,1](31)В нашем случае g(x) это есть функция F (ibx).

Для того, чтобы вычислитькоэффициенты Dj (ω) необходимы значения функций Бесселя. Один из способових вычисления указан в [3]. Вычисления можно производить с использованиемсовременных математических пакетов (например, Maple). Эффективный алгоритм нахождения корней многочленов Лежандра и построения формулы Гаусса(28) приведен в классической работе [60]. Если точность оценки (31) для выбранного n недостаточна, то предлагается поступить следующими способами:1) построить формулу Гаусса с большим числом узлов; 2) применить составнуюквадратурную формулу на основе формулы (30).В первом случае для построения формулы Гаусса наряду с методом работы [60] можно использовать алгоритм, разработанный в [5].43Во втором случае алгоритм состоит в следующем — надо разбить отрезокинтегрирования в (27) на m равных частей, а затем на каждом полученном отрезке применить формулу (30) с фиксированным числом узлов n.

Заметим, чтооценка погрешности вида (31) для каждого отрезка станет в mn раз меньше,а полная погрешность вычисления по составной формуле уменьшится в mn−1раз.Как показали численные эксперименты, достаточно брать n 6 25, а требуемой точности вычисления интегралов вида (27) добиваться за счет увеличенияm.Теперь рассмотрим вычисление интегралаZZ ∞ez t F (z) dz = e−ibte−xt F (−x − bi) dx.L1(32)0Первый способ состоит в применении квадратуры типа Гаусса с весомЛагерра e−x . Но этот путь малоэффективен, так как может потребоваться построение формулы с очень большим числом узлов и с высокой точностью.Пусть на линии L1 выполняется неравенство |F (z)| 6 M.

ПоложимZ ∞Z Xe−xt F (−x − bi) dx ≈e−xt F (−x − bi) dx(33)00и выберем число X так, что для заданного положительного числа ε выполняется неравенствоZ ∞Ze−xt F (−x − bi) dx 6 MX∞e−xt dx =XM −Xte6 ε.tНа отрезке [0, X] подынтегральная функция в (33) не является осциллирующей,поэтому для вычисления правой части в формуле (33) можно применить составную формулу Гаусса (28) с n узлами и разбиением отрезка [0, X] на m равныхчастей. Погрешность Rn (v) формулы Гаусса (28) оценивается неравенством22n+1 (n!)4|Rn (v)| 6max |v (2n) (x)|.32n + 1 ((2n)!) [−1,1]Следовательно, погрешность составной формулы Гаусса уменьшится в m2n−1раз.

Таким образом, за счет выбора параметров ε и m можно добиться вычисления интеграла (32) с требуемой точностью. Интеграл по линии L3 .44Метод, описанный в этой главе, был применён для вычисления функцииползучести (пространственная переменная x опущена)Ztε(t) = 1 + 4Э−0.8 (τ )dτ0по её изображениюε(p) =14+.p p(p0.2 + 1)В Приложении Программа 4 представлена программа для вычисленияфункции–оригинала ε(t) по её изображению с помощью замены линии интегрирования на параболический контур.

Для получения нужной точности вычисления в программе изменялись значения параметра c — сдвиг контура, необходимы для того, чтобы все особые точки изображения лежали внутри контураинтегрирования, а также параметр N — число узлов формулы трапеций. Каквидно из приведённых ниже Таблиц 11, 12, для конкретной вычисляемой функции есть пороговое значение c = 0.1, при котором точность вычислений падает,но её можно улучшить, увеличивая число узлов формулы трапеций N .Результаты вычисления ε(t) при α = −0.8,Таблица 11Ncε(t)ε(t)точное100.11.867651.86696100.31.866961.86696500.11.866961.8669650 0.01 1.869051.86696t = 0.001.45Результаты вычисления ε(t) при α = −0.8,t = 1.Таблица 12Ncε(t)ε(t)точное100.13.116693.11559100.33.115593.11559500.13.115593.1155950 0.06 3.115593.11559Теперь приведём результаты вычисления функции sin(t) с помощью Программы 5 из Приложения.

В представленной программе проводились вычисления также с помощью замены линии интегрирования на параболический контур, но для получения нужной точности изменялся сдвиг изображения на величину r для того, чтобы все его особые точки лежали внутри контура интегрирования. В конце программы для восстановления функции–оригинала происходитдомножение на ert .Результаты вычислений представлены в Таблицах 13, 14.Результаты вычисления sin(t) при t = 2.Таблица 13Nrsin(t)прибл sin(t)точное100 0.650.909780.90929500 0.650.909290.909291000.70.909790.909295000.70.909290.9092946Результаты вычисления sin(t) при t = 200.Таблица 14Nrsin(t)прибл sin(t)точное1000.3−0.87321−0.873295000.3−0.87321−0.87329100 0.31 −0.87266−0.87329500 0.31 −0.87266−0.87329В Приложении Программа 6 представлена программа для вычисленияфункции–оригинала ε(t) по её изображению с помощью замены линии интегрирования на кусочно прямолинейный контур.

Для получения нужной точностивычисления в программе изменялись значения параметра b — расположениеверхнего и нижнего участка контура, подбирался из условия, чтобы все особыеточки изображения лежали внутри контура интегрирования. С этой же цельюв программе задаётся параметр r — сдвиг изображения, а также параметр m1— число узлов составной квадратурной формулы Гаусса и m — число равныхчастей, на которые разбивается отрезок интегрирования [0, X].Результаты вычисления ε(t) приα = −0.8, t = 0.001, b = 0.001, r = 0.24 · 10−5 .Таблица 15m m1ε(t)ε(t)точное30 20 1.423501.8669640 30 1.496951.8669665 45 1.882861.8669647Результаты вычисления ε(t) при α = −0.8, t = 1, b = 0.14, r = 0.003.Таблица 16m m1ε(t)ε(t)точное30 20 3.115673.1155940 30 3.115673.1155965 45 3.115593.11559Как видно из Таблиц 15, 16 использование кусочно прямолинейного контура для вычисления функции–оригинала по её изображению не является оптимальным, ввиду трудностей, связанных с подбором значений параметров b, rдля получения наилучшего результата вычислений.

Заметим, что подбор параметров зависит от значения аргумента вычисляемой функции. Так, в случаеt = 1 удалось найти b, r, для которых результаты слабо зависят от m и m1(Таблица 16), а в случае малого значения аргумента t = 0.001 в Таблице 15проиллюстрирована сильная зависимость результатов вычислений от m и m1 .Таким образом, данный метод можно рекомендовать для вычисления функции–оригинала лишь при значениях аргумента t > 1.482.2.1.Задачи линейной вязкоупругостиПостановка задач линейной вязкоупругости.Одной из областей применения интегрального преобразования Лапласаявляются задачи линейной вязкоупругости. При малых деформациях, если убрать нагрузку, деформации исчезают и физические тела проявляют упругиесвойства.

Так в случае абсолютно упругого тела напряжение, то есть отношениемодуля силы упругости к площади поперечного сечения, подчиняется законуГукаσ = Eε,где E — это модуль упругости, который характеризует упругие свойства материала.Однако, в реальности материалы отклоняются от закона Гука, например,проявляя и вязкоподобные, и упругие свойства.

Так полимеры, даже при комнатных температурах проявляют способность медленно деформироваться вовремени при постоянных напряжениях. С течением времени ползучесть можетпривести к значительным изменениям в напряженно–деформированном состоянии конструкций.Экспериментально ползучесть можно исследовать при растяжении образцов материала. Верхний конец образца закрепляется, а к нижнему — прикладывается нагрузка. Проводится наблюдение за тем, как изменяется длина образцаи строится кривая изменения деформации ε от времени t — кривая ползучести.

При мгновенном приложении нагрузки возникает мгновенная деформацияε0 = σ0 /E0 .С ползучестью материала связана релаксация — явление перераспределения напряжения в теле при постоянной деформации. Пусть, образец "мгновенно" растянули так, что длина стала равной l. Следовательно, была приложенанекоторая сила. После этого образец был закреплен в растянутом состоянии навремя t. Затем его освободили и приложили нагрузку для растяжения образца49до длины l. В процессе экспериментов оказалось, что нагрузка во втором случае потребуется меньше первоначальной. Это связано с явлением релаксацииматериала — уменьшение напряжения при постоянной деформации.В реальных элементах конструкций ползучесть и релаксация материалапроявляются одновременно и взаимосвязанно. Эта связь описывается соотношением Больцмана–Вольтерра, которое является обобщением закона Гука.

Дляпростоты будем говорить об одномерном случае (вязкоупругий стержень).Ползучесть материала описывается соотношением:Zt1 ε(x, t) =σ(x, t) + (E0 − E∞ ) K(t − τ )σ(x, τ ) dτ  ,E0(34)0где K(t) — ядро ползучести, E0 — мгновенный модуль упругости, E∞ — длительный модуль упругости.Релаксация материала описывается соотношением:ZtR(t − τ )ε(x, τ ) dτ,σ(x, t) = E0 ε(x, t) − (E0 − E∞ )(35)0где R(t) — ядро релаксации. Оно является резольвентой ядра ползучести.

Одноядро может быть выражено через другое.Таким образом, соотношение Больцмана–Вольтерра связывает две физические величины — напряжение и деформацию, то есть описывает характеристики материала. Полная же постановка динамической задачи вязкоупругостивключает в себя ещё уравнение движения, начальные и граничные условия.Так, в случае задачи о распространении полубесконечного импульса нагрузкив полубесконечном вязкоупругом стержне уравнения имеют вид [4]:∂ 2u∂σ(x, t) = ρ 2 (x, t),∂x∂t∂uu(x, 0) =(x, 0) = 0,∂tσ(0, t) = σ0 η(t),50где u(x, t) — абсолютное перемещение (ε(x, t) = ∂u(x, t)/∂x), ρ — линейнаяплотность материала стержня, η(t) — единичная функция Хевисайда, возникающая из условия, что к концу одномерного вязкоупругого стержня в моментвремени t = 0 прилагается нагрузка P (t) = σ0 η(t).В качестве ядра релаксации возьмём ядро Ржаницына [4]:aβ β−1R(t) =texp(−at),Γ(β)где a > 0, 0 < β < 1.

Заметим, что при β = 1 ядро Ржаницына переходит вэкспоненциальное ядро.Покажем, как исходная задача о распространении полубесконечного импульса нагрузки в полубесконечном вязкоупругом стержне сводится к задаче вобразах по Лапласу. По определению ε(x, t) = ∂u(x, t)/∂x, тогда из уравнениядвижения для рассматриваемой задачи получим∂ 2σ∂ 2ε(x, t) = ρ 2 (x, t).∂x2∂tПрименяя преобразование Лапласа по t к этому уравнению (далее зависимостьот пространственной переменной x опущена), придём к соотношению000σ x2 (p) = ρ(p2 ε(p) − pε(0) − ε (0)),также из (26) имеемσ(p) = E0 − (E0 − E∞ )R(p) ε(p) =E0 − (E0 − E∞ )aβ(p + a)βε(p).0Учитывая, что ε(0) = ε (0) = 0, получим уравнение, из которого найдём выражение для образа по Лапласу напряжения:ρp2σ x2 (p) −= 0.aβE0 − (E0 − E∞ ) (p+a)β00Решая это дифференциальное уравнение с учётом ограниченности напряженияσ = σ(x, t) при x → +∞, получим выражение для изображения по Лапласу51напряженияσ(p) =σ0expp√!px ρE01/2a)β )(1 − (1 − E∞ /E0 )(a/p +.Поскольку импульс имеет конечную скорость распространения, то в точку сqкоординатой x напряжение приходит не сразу, а за время τx = x/c0 , c0 = Eρ0 ,поэтому положим [4]σ(x, τ + τx ) =0,τ < 0,Ф(x, aτ /ω), τ > 0,где ω = (1 − E∞ /E0 )−1/β , тогда (далее зависимость от пространственной переменной x опущена)ωσ(p) = e−pτx Ф(pω/a),aоткуда ap aФ(p) = expτx σ(ap/ω).ωωПодставив в это выражение σ(p), получим изображение функции Ф(t)αp,Ф(p) = B(p) exp −ωϕ(p)гдеB(p) = σ0 /p,ββ/2ϕ(p) = (p + ω) − 1 + (p + ω)pα = ax ρ/E0 = aτx .q(p + ω)β − 1,Из вида изображения функции Ф(t) следуют предельные соотношения:1) σ(x, t) → σ0 , при t → +∞;2) при 0 < β < 1 σ(x, τx + 0) = 0, в случае же β = 1 σ(x, τx + 0) =axσ0 exp − 2ωc0 .