Поскольку в общем случае поиск точных значений корней на компьютере невозможен, вычисление корней функции осуществляется с заданной точностью.
Подставив
в общий вид функции
параметры для слоя и покрытия, с учетом того, что
- параметры покрытия, а
- параметры слоя, был получен общий вид уравнения,
нули и полюса которого необходимо найти.
Отдельно
выделяются вещественная и мнимая части функции с помощью функций Real
и Aimag. Нули и той и другой частей ищутся на интервале
, для каждого конкретного значения
. Причем
изменяется
от 0 до 10 с шагом, равным 0.05, в то время как для
он равен 0.01. Результаты работы программы выводятся
в отдельные файлы.
Для
поиска нулей использовалась стандартная подпрограмма библиотеки IMSL для поиска
вещественных корней вещественной функции ZBREN [9].
Подпрограмма ZBREN применяет метод бисекций, линейную интерполяцию
(метод секущих) и обратную квадратичную интерполяцию. На каждом шаге
принимается решение, какой из трех методов будет использован для вычисления
следующего приближения. Также IMSL использует для определения заданного числа
вещественных корней метод Мюллера. Последний, вдобавок, применяется и при
вычислении комплексных корней. (Данный метод является трехшаговым и основан
примерно на тех же идеях, что и метод обратной квадратичной интерполяции а
также метод секущих. Геометрически, на каждом шаге поиска очередного
приближения корня используется три точки и строится парабола, проходящая через
эти точки. В качестве следующего приближения берётся точка пересечения параболы
и оси
).
Суть
алгоритма заключается в исследовании отрезка
, на
котором
. Если найден такой промежуток, то запускается
стандартная подпрограмма ZBREN, которая имеет вызов:
Call ZBREN(f, errabs,errel,a,b,maxfn).
Параметры подпрограммы ZBREN:
- Пользвательская
функция:
;
- Входные данные: errabs,errel;
- Входные/выходные данные: a,b,maxfn.
Здесь
errabs - первый критерий завершения вычислений. Корень b
принимается, если
. Причем, можно задать
. Далее, errel - относительная ошибка
(второй критерий завершения вычислений). Корень принимается, если относительная
разница между приближениями, найденными в двух последовательных аппроксимациях,
меньше errel.
Параметры а и b - начало и конец отрезка, на котором выполняется поиск вещественного корня. На выходе они отличны от начального значения, причем в b содержится приближенное значение корня функции f.
Параметр maxfn - на входе задает максимально допустимое число вычислений функции f, на выходе содержит число реальных вызовов функции f.
Значения корней для вещественной и мнимой частей сравниваются и, если они совпадают, то используются в дальнейшем для построения графиков.
Для
анализа полученных данных в среде Maple15 были построены графики
вещественных нулей и полюсов в зависимости от частоты
.
Рисунок
2 - Кривые нулей и полюсов функции
На
рисунке 2 приведены кривые нулей и полюсов функции
(красным цветом обозначены нули элемента, черным -
полюса) при следующих параметрах задачи:
;
;
;
;
;
;
;
.
Также были построены графики, отражающие
расположение нулей и полюсов заданной функции, для различных наборов таких
параметров.
Рисунок 3 - Кривые нулей и полюсов функции для значений
;
;
;
;
;
;
;
Рисунок 4 - Кривые нулей и полюсов функции для значений
;
;
;
;
;
;
;
Рисунок 5 - Кривые нулей и полюсов функции для значений
;
;
;
;
;
;
;
Рисунок
6 - Кривые нулей и полюсов функции для значений
;
;
;
;
;
;
;
Рисунок
7 - Кривые нулей и полюсов функции для значений
;
;
;
;
;
;
;
Из рисунков 2 - 7 видно, что наблюдается чередование вещественных нулей и полюсов, с увеличением частоты число вещественных полюсов растет.
При построении графиков варьировались параметры, характеризующие толщину полосы и покрытия, а также их плотность.
Выразим для рассматриваемых вариантов сред с покрытием интегральные характеристики напряжений и перемещений на стыке покрытия и слоя через известные величины.
Амплитуды напряжений и перемещений могут быть выражены через их
Фурье-образы следующим образом:
,
,
где Г - контур, совпадающий с вещественной осью везде, за исключением отрезка конечной длины, содержащего вещественные полюсы подынтегральной функции. Так как в нашем случае вещественные полюсы однократные, то согласно принципу предельного поглощения, контур Г обходит отрицательные полюсы сверху, положительные - снизу.
Для
вычисления перемещений была создана программа на языке Fortran, которая
производит обратное преобразование Фурье. Использовалась функция DQDAGS из библиотеки
IMSL [9], вычисляющая значение интеграла от функции,
имеющей конечное число особенностей на заданном промежутке. Она использует
адаптивную схему, уменьшающую абсолютную ошибку. Подпрограмма делит отрезок
на подынтервалы и использует 21-точеченое правило
Гаусса-Кронрода для оценки интеграла на каждом подынтервале.
DQDAGS имеет вызов:
Call DQDAGS(F,a,b,Errabs,Errel,Res,Errest).
Входными параметрами DQDAGS являются:
- F - пользовательская функция, интеграл от которой должен быть вычислен;
- а - нижняя граница интегрирования.
- b - верхняя граница интегрирования;
- Errabs - желаемая абсолютная ошибка;
- Errel - желаемая относительная ошибка.
Выходные параметры:
- Res - значение вычисленного интеграла;
- Errest - оценка величины абсолютной ошибки.
Так как подынтегральная функция является комплексной, то вычисления ведутся для вещественной и мнимой части.
На
рисунках 8 - 16 представлены графики вещественной части вертикальных смещений
пластины для различных значений частоты
и для
следующих параметров задачи:
;
;
;
;
;
;
;
.
Рисунок
8 - Вещественная часть вертикальных смещений u
при
=1,9
Рисунок
9 - Вещественная часть вертикальных смещений u
при
=2,4
Рисунок
10 - Вещественная часть вертикальных смещений u
при
=3,4
Рисунок
11 - Вещественная часть вертикальных смещений u
при
=4,4
Рисунок
12 - Вещественная часть вертикальных смещений u
при
=5,4
Рисунок
13 - Вещественная часть вертикальных смещений u
при
=6,6
Рисунок
14 - Вещественная часть вертикальных смещений u
при
=7,3
Для
вычисления напряжений также была написана программа на языке Fortran,
которая производит обратное преобразование Фурье. Принцип ее работы аналогичен
предыдущему. Ниже приведены графики, построенные по результатам работы
программы для различных значений частоты
.
Значения вычислялись при следующих значениях параметров:
;
;
;
;
;
;
;
.
Рисунок
17 - Вещественная часть амплитуды напряжений q
при ????=1,4
Рисунок
18 - Вещественная часть амплитуды напряжений q
при
=2,4
Рисунок
19 - Вещественная часть амплитуды напряжений q
при
=3,6
Рисунок
20 - Вещественная часть амплитуды напряжений q
при
=4,2
Рисунок
21 - Вещественная часть амплитуды напряжений q
при
=5,3
Рисунок
22 - Вещественная часть амплитуды напряжений q
при
=6,6
Рисунок
23 - Вещественная часть амплитуды напряжений q
при
=7,6
Рисунок
24 - Вещественная часть амплитуды напряжений q
при
=8,3