| МИНИСТЕРСТВО ОБРАЗОВАНИЯ РОССИЙСКОЙ ФЕДЕРАЦИИ
КУБАНСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
Кафедра математического моделирования
КУРСОВАЯ РАБОТА
ИССЛЕДОВАНИЕ ДИНАМИКИ УПРУГОГО ПРОСТРАНСТВА, СОДЕРЖАЩЕГО СИСТЕМУ ПЛОСКИХ ВКЛЮЧЕНИЙ
Работу выполнила Вербенко Н.П.
Группа 43 Факультет прикладной математики Специальность 0102
Руководитель работы к.ф.- м.н., доцент Рубцов С.Е.
Краснодар 2002
РЕФЕРАТ
Курсовая работа посвящена изучению методов решения краевых задач динамической теории упругости для неограниченных выпуклых областей с неплоской границей, частных случаев включений с плоской границей и ряда численных методов для применения к решению задач для сред с неоднородностями. А также проведению вычислений для плоских включений
Курсовая работа содержит 48 страниц, 4 рисунка и 3 приложения.
СОДЕРЖАНИЕ
Введение 4
1. Задачи теории упругости 5
1.1. Вывод интегральных уравнений для полупространства 9
1.2. Полупространство с плоской границей 14
2. Постановка задачи для пространства с внутренними плоскими включениями 16
3. Решение задачи для пространства
3.1 Случай одного включения 19
3.2 Случай двух включений 24
3.3 Случай трех включений 27
4. Поиск вещественных полюсов
4.1 Случай одного включения 31
4.2 Случай трех включений 32
Заключение 33
Список литературы 34
Приложение 1 35
Приложение 2 39
Приложение 3 43
ВВЕДЕНИЕ
Динамические эффекты в различных средах в связи со своей теоретической и практической значимостью для различных областей деятельности человека стали на сегодняшний день одним из основных объектов изучения. Их исследование способствует решению многочисленных задач сейсмологии, физики Земли, машиностроения и техники в целом.
В данной работе описаны постановки смешанных задач теории упругости для неоднородностей в упругих изотропных средах и решена задача для пространства для случая одного и двух включений. В качестве неоднородностей рассматриваются включения, для которых соотношения между напряжениями и перемещениями на границах сводятся к системам интегральных уравнений, которые в свою очередь преобразуются к матричному представлению. Для исследования распределения нулей определителей в зависимости от различных условий на параметры включения написана программа на языке С+.
+. Приведены некоторые из полученных расчетов.
1. ЗАДАЧИ ТЕОРИИ УПРУГОСТИ
Колебания упругой среды (в перемещениях) описываются уравнениями Ляме:
 , (1.1)
где λ, μ
– константы Ляме, положительность которых обеспечивает обратимость закона Гука, ρ
– плотность среды.
Вектор  (x1
,
x2
,
x3
,
t
)=(u(
x1
,
x2
,
x3
,
t),
v(
x1
,
x2
,
x3
,
t),
w(
x1
,
x2
,
x3
,
t))
T
– вектор перемещений.
Рассматриваемая упругая среда может быть пространством, полупространством (в том числе стратифицированным), слоем или пакетом слоев.
Механическое состояние упругого тела, занимающего в начальном состоянии известный объем V
с ограничивающей поверхностью S
, характеризуется компонентами тензора деформации ε
ij
и тензора напряжений σ
ij
.
Перемещения в точках тела, под действием заданной системы поверхностных и объемных сил, описываются вектором перемещений
 
и представляют собой непрерывные и однозначные функции координат и времени. Механическое состояние упругого тела характеризуется также вектором напряжений τ ,
возникающих в упругом теле на некоторой элементарной площадке с нормалью  . Вектор τ
выражается через компоненты тензора напряжений  .
Вектор напряжений τ
можно выразить через перемещения  =
T
,
где Т
– линейный дифференциальный оператор напряжений. В изотропном случае
 .
В случае установившихся колебаний зависимость всех характеристик задачи от времени t
описывается множителем  , а уравнения линейной упругости, описывающие установившийся закон колебаний, имеют вид:
 , (1.2)
Для постановки задач теории упругости должны быть заданы, индивидуальные для конкретных задач, граничные условия, а в случае динамических задач – еще и начальные условия.
Интегральное представление решения задачи строится на основе интегральных соотношений между напряжениями и перемещениями, существующих для упругих сред. При этом используется формула Бетти в виде
 (1.3)
здесь 
В этом соотношении вектор–функции u
, v–
произвольные дважды непрерывно дифференцируемые. Область W
может быть неограниченной с гладкой границей.
Вводится система следующих векторов
 (1.4)
Эти векторы формируют матрицу
 (1.5)
В соотношении (1.2) вместо произвольного вектора u
вносится решение дифференциального уравнения  (1.6)
, а вместо вектора v
–матрицу V
. Тогда соотношение (1.2) порождает не скалярное, а векторное выражение.
Выражение  имеет вид
 (1.7)
Матрицу справа обозначим как В
, тогда векторы  – формируемые строками этой матрицы. Можно записать
 (1.8)
Вычисляются векторы  . В результате вычислений будут получены формулы
 (1.9)
здесь для m=1, 2, 3
берется соответственно по порядку n=2 ,
p=3;
n=1,
p=3;
n=1,
p=2
; 
Внося в эти соотношения вместо u
значения векторов  , получаем
 (1.10)
Вектор напряжений, действующих на границе S,
будем обозначать
 (1.11)
где Tk
–
проекции вектора напряжений на направления осей xk
соответственно.
Вычисляя вектор Tu
на границе S
с учетом (1.9) и соотношения
получаем выражение для вектора Tu
:
Tu=
T0
(1.12)
Внося соотношения (1.9) , (1.10) , (1.12) , в (1.2) для каждого вектора, и вектора решения уравнения (1.2) приходим к соотношению
 (1.13)
Где I
–единичная матрица третьего порядка.
Если разрешить соотношения (1.13) относительно тройного интеграла в левой части, то получим преобразование Фурье искомого решения краевой задачи для дифференциального уравнения (1.6). Оно будет зависеть как от значения вектора перемещений на границе области, так и от вектора напряжений. Так как в любой области на границе S
может быть задан согласно постановке задачи теории упругости лишь один из этих векторов, то между компонентами вектора перемещений и напряжений на границе должны существовать определенные зависимости.
1.1. ВЫВОД ИНТЕГРАЛЬНЫХ УРАВНЕНИЙ ДЛЯ ПОЛУПРОСТРАНСТВА
Для вывода интегральных уравнений для полупространства с выпуклой границей рассмотрим следующую матрицу
 (1.14)
Ее определитель D
после вычислений принимает вид
 (1.15)
Вычисляем обратную к (1.14) матрицу
 (1.16)
В дальнейшем элементы матрицы (1.16) будут обозначаться ckm
.
Таким образом, имеем
 (1.17)
Теперь подействуем слева на соотношение (1.13) матрицей (1.16), после вычисления найдем
 (1.18)
p=1,2,3
Для вывода интегральных уравнений, связывающих на границе перемещения um
и напряжения Tk
нужно изучить аналитические свойства функции  .
Учитывая геометрию тела, заключаем, что она должна быть аналитически продолжаемой по параметру a3
в верхнюю полуплоскость. Это вытекает из свойства области, описываемой неравенством
Для выяснения условия, обеспечивающего аналитическую продолжимость решения, изучим распределение особенностей в представлении функции Up
. Все особенности этих функций описываются уравнением
 (1.19)
Находим корни этого уравнения
  (1.20)
Но  . Таким образом, особенности функций описываются корневыми множествами функций трех комплексных переменных вида
 (1.21)
Эти корневые множества представляются аналитическими подмногообразиями в C3
комплексной размерности, равной двум. Разрешая уравнение относительно a3
, получаем представление корневых множеств вида
Здесь у радикала выбраны такие ветви, что
 при  (1.22)
В этом случае при  имеем
 (1.23)
Поскольку область содержит положительную полуось, то функции Up
не должны иметь особенностей по параметру a3
в верхней полуплоскости, т.е. должны быть аналитически продолжимыми в область  .
Следовательно, правая часть в выражении (1.18) должна быть ограниченной при  и  .
Для построения вытекающих из этого условия соотношений положим в (1.18):
 (1.24)
Тогда для ограниченности Up
достаточно потребовать обращения в нуль при и  выражений (1.18), в которых вместо  взято  .
Выполнив эти подстановки, приходим к выражениям вида:
 (1.25)
 (1.26)
Здесь lk
–функция параметров.
Соотношения (1.25) , (1.26) являются своеобразной формой записи интегральных уравнений, связывающих заданные на границе S
тела напряжения Tk
и перемещения um
.
Так, если рассматривается краевая задача теории упругости I рода, т.е. при заданных на границе S
напряжениях, то правые части соотношений (1.25) , (1.26) - известные выражения. И из уравнений необходимо найти перемещения um
, стоящие в левых частях.
Если же рассматривается краевая задача теории упругости II рода, т.е. при заданных на границе S.
перемещениях, то, наоборот, левые части соотношений (1.25) , (1.26) известны, а определению подлежат стоящие в правых частях функции Т
k
.
В случае смешанной задачи теории упругости, когда на одном множестве S1
поверхности S
заданы напряжения Т
k
, а на другом S2
-
перемещения и
m
,
соотношения (1.25), (1.26) представляют собой систему интегральных уравнений, неизвестные которых расположены и в левой, и в правой частях в зависимости от того, по какому множеству осуществляется интегрирование. Важно, что на одном и том же множестве не могут быть одновременно либо только неизвестные и справа и слева, либо только известные функции.
Допустим, что удалось решить интегральные уравнения (1.25), (1.26) для перечисленных выше краевых задач теории упругости. Тогда известными функциями на всей поверхности S
оказываются как перемещения и
m
,
так и напряжения Tk
.
. Внесем их в (1.18).
Геометрия задачи (полупространство) позволяет применить преобразование Фурье. Вообще преобразования Фурье при решении краевых задач теории упругости наиболее эффективно в тех случаях, когда упругое тело занимает объем, содержащий бесконечно удаленные точки, и все границы тела параллельны тем координатным осям, по которым берется преобразование. Свойства тела вдоль этих направлений также должны быть постоянными, но могут меняться в перпендикулярных направлениях. Преобразование Фурье не применимо по тем координатам, вдоль которых свойства среды непрерывно меняются.
В результате выражение для up
(
x1
,
x2
,
x3
)
после применения формул обращения Фурье можно представить в виде:
 (1.27)
Интеграл по параметру a3
.
вычисляется по теории вычетов.
Из соотношения (1.11) , (1.12) видно, что решение содержит волны
уходящие от поверхности S
на бесконечность.
1.2. ПОЛУПРОСТРАНСТВО С ПЛОСКОЙ ГРАНИЦЕЙ
Рассмотрим частный случай исследуемой задачи, когда поверхность  вырождается в плоскость, совпадающую с координатной плоскостью x3
=0
. В этом случае внешняя нормаль к поверхности остается неизменной в любой точке и имеет компоненты
 (1.28)
Кроме того, надо принять x3
=0 на S
, а также равенства
 (1.29)
Внесем значения указанных параметров в интегральные уравнения (1.11), (1.12). Тогда в подынтегральных выражениях члены
оказываются постоянными, не зависящими от параметров  и их можно вынести за знак интегралов. Введем следующие обозначения:
 (1.30)
Интегральные уравнения (1.11), (1.12) можно переписать в виде
 (1.31)
Эти соотношения полностью совпадают с уравнениями, получающимися при удовлетворении граничных условий в краевой задаче об установившихся колебаниях упругого полупространства. Если на границе задается вектор напряжений То
, то для получения перемещений на границе левые части соотношения (1.31) необходимо разрешить относительно компонент вектора U0
.
В матричном виде систему (1.31) можно представить в форме
 (1.32)
2. ПОСТАНОВКА ЗАДАЧИ ДЛЯ ПРОСТРАНСТВА
С ВНУТРЕННИМИ ПЛОСКИМИ ВКЛЮЧЕНИЯМИ
Упругая изотропная среда содержит систему произвольного количества – N
плоских включений, расположенных параллельно на высотах h1
,…,hN
соответственно (h1
<…< hN
). Включения занимают области Ωl
с границами Sl
, внешние нормали к границам  , l=1,…,N
.
Колебания упругой среды описываются уравнениями Ляме (1.1). Вектор перемещения точек среды в данном случае  (x1
, x2
, x3
, t)
.
Дополнительные условия, налагаемые на рассматриваемые характеристики на границах заданных областей, приводят к краевым задачам. Упругая среда может быть пространством, полупространством или слоем. При этом на верхней границе могут быть заданы компоненты вектора перемещений
или вектора напряжений  .
Для установившихся колебаний (
x1
,
x2
,
x3
,
t)
= (
u1
,
u2
,
u3
)е-
iωt
, уравнение Ляме (1.1) относительно комплексной амплитуды u
примет вид:
 . (2.2)
Для массива включений, на берегах неоднородностей заданы перемещения:
 (2.3)
При этом напряжения среды в области включений терпят разрыв:
 (2.4)
В покомпонентной записи уравнение (1.1) имеет вид
 (2.5)
Для упругого пространства, содержащего плоские неоднородности, в качестве граничных условий берутся перемещения берегов включений:
 , l =
1 ,…, N
(2.6)
и условия убывания на бесконечности
 , (2.7)
дополненные условиями излучения.
В качестве условий излучения выбраны следующие принципы.
1. Принцип Зоммерфельда: в решении удерживаются составляющие, описывающие волны, уходящие от источника в бесконечность, и отбрасываются те, скорость которых направлена к источнику.
2. Принцип предельного поглощения: в качестве решения задачи для идеально упругой среды берется равномерный по всем параметрам предел решения соответствующей задачи для вязкоупругой среды (среды с поглощением) при стремлении вязкости к нулю.
3. РЕШЕНИЕ ЗАДАЧИ ДЛЯ ПРОСТРАНСТВА
3.1.СЛУЧАЙ ОДНОГО ВКЛЮЧЕНИЯ
Разделим пространство по плоскому включению на два полупространства. Для каждого из них справедливо соотношение (1.32). Введем обозначения:
Тв
,
Uв
–
напряжение и перемещение сверху включения,
Тн
,
Uн
–
напряжение и перемещение снизу включения.
Тогда задача для одного включения сведется к системе функциональных соотношений:
 (3.1)
в предположении, что
Uв
=
Uн
=
U
(3.2)
Учитывая соотношения из раздела 2 , L+
,
L–
,
D+
,
D
–
принимают вид:
 (3.3)
Выяснив явный вид матриц в соотношениях (3.1), получим напряжение на границе плоскости:
Tв
=(
D+
)-1
L+
Uв
(3.4)
Tн
=(
D–
)-1
L–
Uн
Учитывая, что для данного случая с одним включением имеем следующие представления корневых множеств:
 (3.5)
а так же c
1
= , c
2
= .
 Введем обозначение
 (3.6)
Решим уравнение (3.4). Найдем матрицы, обратные к D+
и D
–
.Для этого:
а) Вычисляем определители матриц
 (3.7)
где  (3.8)
б) Транспонируем матрицы D+
и D
–
.
 (3.9)
в) В соответствии с теорией матриц [2], определяем дополнительные миноры транспонированных матриц и выписываем обратные матрицы:
 (3.10)
 (3.11)
Умножив полученные обратные матрицы (3.10), (3.11) на исходные матрицы L+
и L–
соответственно, найдем
 (3.12)
где  ,
 ,
 ,
 ,
 ,
 ,
 ,
 ,
 ,
и
 , (3.13)
где
 ,
 ,
 ,
 ,
 ,
 ,
 ,
 .
Итак, если на границе включения задается вектор перемещений U, то скачок вектора напряжений на границе имеет следующий вид:
T+
– T–
=i
m
-
1
k
-
1
MU
, (3.14)
где матрица M – разность произведений матриц (
D+
)-1
L+
и (
D–
)-1
L
–
, выглядит следующим образом:
 (3.15)
Для получения вектора перемещений из (3.14)
 ,
при заданном скачке напряжений Т=Т+
– Т
–
необходимо найти обратную матрицу к матрице М.
Введем обозначение:
Определитель матрицы М равен:
 ,
где
 (3.16)
 (3.17)
Обратная к М
матрица равна:
 (3.18)
3.2. СЛУЧАЙ ДВУХ ВКЛЮЧЕНИЙ.
Пусть упругое изотропное пространство ослаблено двумя плоскостями, расположенными на высотах h1
и h2
,моделирующими включения бесконечно малой толщины. Сохраняя все остальные предположения и обозначения из пункта 3.1, решение данного частного случая проведем аналогичным образом.
Матричный вид системы, удовлетворяющей граничным условиям в краевой задаче об установившихся колебаниях упругого пространства, будет следующим:
 (3.21)
в предположении что  . Матрицы  имеют вид:
 (3.22)
 ,
где m=1,2.
Таким образом, если на границах включений задаются вектора напряжений  , то для получения перемещений на границе систему (3.21) необходимо разрешить относительно компонент U1
и U2
.
Введем обозначения:
 (3.23)
После несложных действий над матричной системой (3.21) с учетом правил умножения матриц приходим к системе вида:
 (3.24)
С учетом (3.23) в итоге получаем
 (3.25)
где  дается соотношением (3.15),
 (3.26)
Они имеют следующий вид:
 (3.27)
(3.28)
где x=
iα31
(
h2
-
h1
)
и y=
iα32
(
h2
-
h1
)
.
3.3 СЛУЧАЙ ТРЕХ ВКЛЮЧЕНИЙ
Пусть упругое изотропное пространство ослаблено тремя плоскими включениями, расположенными на высотах h1
,
h2
,
h3
.Сохраняя все предположения и обозначения из (3.1), решение данного частного случая проведем аналогичным образом.
Матричный вид системы будет следующим:
 (3.31)
в предположении, что  , m=1,2,3.
Матрицы  имеют вид (3.22), где m=1,2,3.
Т.о., если рассматривается краевая задача теории упругости I рода, т.е. при заданных на границе S напряжениях, то для получения перемещений на границе систему (3.31) необходимо разрешить относительно компонент U1
,
U2
,
U3
.
Обозначим  .
После несложных преобразований с учетом правил умножения матриц и введенных обозначений приходим к системе вида:
 , (3.32)
В общем виде система матричных уравнений (3.32) имеет вид:
 (3.33)
где матрицы  ,   имеют вид (3.18), (3.27) и (3.28) соответственно.
Вычисление определителей:
Вычисление определителей матриц можно провести на основе метода Гаусса. Этот метод является наиболее простым и естественным для решения систем уравнений и основан на последовательном исключении неизвестных. Метод Гаусса имеет много различных вычислительных схем. Остановимся на описании соответствующей схеме единственного деления.
Пусть
 (3.34)
и пусть a11
 .
Вынесем элемент из первой строки. Тогда, применяя следующие обозначения  , получим:
Затем из каждой строки отнимаем первую строку, умноженную соответственно на первый элемент этой строки. Очевидно, получим следующее:
Далее с оставшимся определителем (
n-1)
порядка поступаем совершенно таким же образом, если только  .
Продолжая процесс, получим, что искомый определитель равен произведению ведущих элементов:
Если на каком-либо шаге окажется, что  или  близко к нулю (что влечет за собой уменьшение точности вычислений), можно предварительно изменить порядок строк и столбцов определителя так, чтобы в левом верхнем углу оказался неисчезающий элемент. Наилучший в смысле надежности результат получится, если на каждом шагу процесса переводить в левый верхний угол наибольший по абсолютной величине элемент матрицы, для которой вычисляется определитель – схема единственного деления по главным элементам.
Число умножений и делений, нужных для вычисления определителя nго
порядка, равно (
n-1)(
n2
+
n+3)/3
.
Схема единственного деления для вычисления определителей может применяться не только с исключением по строкам, но и по столбцам. Это равносильно вычислению определителя транспонированной матрицы по строкам. Получение нулей при вычислении определителя можно осуществить посредством линейного комбинирования строк матрицами вращения. Это требует значительно большего числа операций, но дает значительно более стабильный процесс. Для вычисления определителя можно так же разбить матрицу на клетки с квадратными диагональными клетками и затем “получить нули”. Разбиение матрицы на клетки, вообще говоря, не изменяет числа элементарных операций, необходимых для вычисления, но можно получить значительный выигрыш в объеме работы. В конце концов, искомый определитель окажется произведением определителей “ведущих” клеток.
4. ПОИСК ВЕЩЕСТВЕННЫХ ПОЛЮСОВ
4.1 Случай одного включения
Для получения решения исходной задачи необходимо выполнить обратное преобразование Фурье.
Одним из способов вычисления обратных интегралов Фурье является метод, основанный на использовании теории вычетов. Применение этого метода требует исследования функции
и
на наличие нулей для различных значений частоты w
.
Нули функции   можно найти аналитически, а для исследования   необходимо проводить численные расчеты.
В результате несложных преобразований для найдены нули, которые равны  .
Для численного нахождения нулей функции  написана программа на языке С+ +, алгоритм которой заключается в следующем:
Для различных значений частоты w, изменяющихся в пределах от 1 до 10 включительно с шагом 0.1, берутся различные b
и методом дихотомии находятся такие значения b
, при которых функция ∆2
обращается в нуль.
Исходя из равенств (3.5), берем во внимание только те значения β
, которые удовлетворяют следующим условиям:
b<
c1
c2
<
b<2
c2
,
т.к. на интервале (
c1
,
c2
)
нули функции не будут вещественными.
Численные расчеты по данной программе показали то, что функция ∆2
не имеет вещественных нулей в указанном интервале.
4.2 Случай трех включений
Аналогично случаю одного включения, для решения краевой задачи об установившихся колебаниях упругого пространства необходимо вычислить и исследовать на наличие вещественных нулей определитель блочной матрицы коэффициентов системы (3.33):
Для проведения анализа поведения функции определителя матрицы С
была написана программа, в которой для различных значений a1
и a2
, изменяющихся в диапазоне от 0
до  , и различных высот  включений находится определитель и строится график зависимости  от a1
и a2
.
ЗАКЛЮЧЕНИЕ
В данной работе исследовалась упругая задача для сред с включениями, напряженно–деформированное состояние которых определялось решением систем уравнений как относительно перемещений, так и относительно напряжений.
Было получено решение в образах Фурье задачи для пространства с одним, двумя и тремя включениями.
Для задач с одним и тремя включениями были найдены нули определителей частью аналитически, частью – численно – с помощью специально написанных программ на языке С+ +. Случай с тремя включениями был исследован для различных высот плоских включений относительно друг друга.
В дальнейшем предполагается выполнение обратного преобразования Фурье, рассмотрение случая многих включений, а также исследование задачи с включениями для полупространства и слоя.
СПИСОК ЛИТЕРАТУРЫ
1. Бабешко В.А Динамика неоднородных линейно-упругих сред/
В.А. Бабешко, Е.В. Глушков, Ж.Ф. Зинченко.– М. Наука. 1989.
2. Гантмахер Ф.Р. Теория матриц – М. Наука. 1988.
3. Фаддеев Д.К. Вычислительные методы линейной алгебры/
Д.К. Фаддеев, В.Н. Фаддеева – М. 1960.
4. Демидович Б.П. Основы вычислительной математики/
Б.П. Демидович, И.А. Марон.– М.Наука.1970.
 ПРИЛОЖЕНИЕ 1
 
  4
ПРИЛОЖЕНИЕ 2
#include <iostream.h>
#include <math.h>
double alpha(double x,double b)
double talpha=sqrt(x*x-b*b);
return talpha;
double alpha1(double d)
double talpha=sqrt(fabs(d));
return talpha;
double s(double x,double b)
double ss=0.5*x-b*b;
return ss;
double f(double b,double a31,double a32,double s)
double fff=a32*a32*b*b*(a31*a31+s)+a31*a32*(b*b*b*b+s*a32*a32);
return fff;
double ff(double b,double a12,double a1,double a2,double s)
double kk=a2*b*b*(a1+s)+a12*(b*b*b*b+s*a2);
return kk;
double compl(double x,double b)
double g=x*x-b*b;
return g;
void koren(double b,double w)
cout<<"При w="<<w<<"есть корень, равный "<<b<<endl;
void main()
{
double p,tt,t,k,r=0.0,bb=0,e=0.001,amega=10.0,betta=0.0,c1=5.6,c2=3.5,c,x,x1=0,x2=0,s1=0;
double h_betta=0.1,h_amega=0.1,func1,func2;
double f(double,double,double,double);
double compl(double,double);
double ff(double,double,double,double,double);
int v=0;
amega=h_amega;
do {
amega+=h_amega;
betta=h_betta;
x1=amega/c1;
x2=amega/c2;
do {
c=betta;
do {
betta+=h_betta;
if (betta>x1) betta=x2+h_betta;
else
{
func1=f(betta, alpha(x1,betta), alpha(x2,betta), s(x2,betta));
if (func1==0)
{koren(betta,amega);v=1;}
}
bb=betta+h_betta;
if (bb>x1)
bb=x2+h_betta;
else
{
func2=f(bb, alpha(x1,bb), alpha(x2,bb), s(x2,bb));
if (func2==0)
{koren(bb,amega);v=1;}
}
if (func1*func2<0)
{x=(betta+bb)/2;
r=f(x,alpha(x1,x),alpha(x2,x),s(x2,x));
if (r*func1<0)
bb=x;betta=bb;
if (r*func2<0) betta=x;
if (r*func1==0)
{koren(betta,amega);v=1;}
}
if (betta>x2&&betta<2*x2)
{
t=compl(x1,betta);
tt=compl(x2,betta);
func1=ff(betta, alpha1(t*tt), alpha1(t*t), alpha1(tt*tt), s(x2,betta));
if (func1==0)
{koren(betta,amega);v=1;}
bb+=h_betta;
if (bb<2*x2)
{
t=compl(x1,bb);
tt=compl(x2,bb);
func2=ff(bb, alpha1(t*tt), alpha1(t*t), alpha1(tt*tt), s(x2,bb););
if (func2==0) {koren(bb,amega);v=1;}
if (func1*func2<0)
{
x=(betta+bb)/2;
r=f(x,alpha1(compl(x1,x)),alpha1(compl(x2,x)),s(x2,x));
if (r*func1<0) bb=x;betta=bb;
if (r*func2<0) betta=x;
if (r*func1==0) {koren(x,amega);v=1;}
}
while (r>e|v==0|betta<2*x2);
betta=c;
}
while (betta<2*x1);
}
while (amega<10);
if (v==0) cout<<"в данном уравнении нет вещественных корней!"<<endl;
}
ПРИЛОЖЕНИЕ 3
#include <fstream.h>
#include <iostream.h>
#include <math.h>
#include <complex.h>
#include <io.h>
#include <process.h>
double alpha31(double x,double a1,double a2)
{ double talpha;
talpha=sqrt(fabs(x*x-a1*a1-a2*a2));
return talpha;}
double s(double x,double a1,double a2)
{ double ss=0.5*x-a1*a1-a2*a2;
return ss;}
double cabs(double x,double y)
{ double g=sqrt(x*x+y*y);
return g;}
complex det(double a1,double a2, double a31,double a32,double s, double \
h1,double h2,double h3)
{ complex det1[9],det2;
complex ii(0,1);
complex x12=exp(ii*ii*a31*(h2-h1));
complex x13=exp(ii*ii*a31*(h3-h1)),x21=exp(ii*ii*a31*(h2-h1)),\
x32=exp(ii*ii*a31*(h2-h3));
complex x23=exp(ii*ii*a31*(h3-h2)),x31=exp(ii*ii*a31*(h1-h3));
complex y12=exp(ii*ii*a32*(h2-h1)),y13=exp(ii*ii*a32*(h3-h1)),\
y21=exp(ii*ii*a32*(h2-h1)),y32=exp(ii*ii*a32*(h2-h3));
complex y23=exp(ii*ii*a32*(h3-h2)),y31=exp(ii*ii*a32*(h1-h3));
complex A[9][9];
complex delta1=s*a31+a32*(a1*a1+a2*a2);
complex delta2=a32*a32*(a1*a1+a2*a2)*(a31*a31+s)+a31*a32*(a1*a1+\
a2*a2+s*a32*a32);
complex M[9][9]={{(a31*(2*a2*a2+a32*a32)+a32*(a1*a1-a2*a2))/delta2,2*a1*a2*\
(a31-a32)/delta2,0,(a31*(2*a2*a2+a32*a32)+a32*(a1*a1-a2*a2))/delta2,2*a1*a2*\
(a31-a32)/delta2,0,(a31*(2*a2*a2+a32*a32)+a32*(a1*a1-a2*a2))/delta2,2*a1*a2*\
(a31-a32)/delta2,0,},\
{a2*(2*a1*a1-a32*a32)/(a1*delta1),2*(a1*a1*a31+a32*(a2*a2+s))/delta2,0,\
a2*(2*a1*a1-a32*a32)/(a1*delta1),2*(a1*a1*a31+a32*(a2*a2+s))/delta2,0,\
a2*(2*a1*a1-a32*a32)/(a1*delta1),2*(a1*a1*a31+a32*(a2*a2+s))/delta2,0},\
{0,0,1/delta1,0,0,1/delta1,0,0,1/delta1},\
{(a31*(2*a2*a2+a32*a32)+a32*(a1*a1-a2*a2))/delta2,2*a1*a2*(a31-a32)/delta2,0,\
(a31*(2*a2*a2+a32*a32)+a32*(a1*a1-a2*a2))/delta2,2*a1*a2*(a31-a32)/delta2,0,\
(a31*(2*a2*a2+a32*a32)+a32*(a1*a1-a2*a2))/delta2,2*a1*a2*(a31-a32)/delta2,0},\
{a2*(2*a1*a1-a32*a32)/(a1*delta1),2*(a1*a1*a31+a32*(a2*a2+s))/delta2,0},\
{0,0,1/delta1,0,0,1/delta1,0,0,1/delta1},\
{(a31*(2*a2*a2+a32*a32)+a32*(a1*a1-a2*a2))/delta2,2*a1*a2*(a31-a32)/delta2,0,\
(a31*(2*a2*a2+a32*a32)+a32*(a1*a1-a2*a2))/delta2,2*a1*a2*(a31-a32)/delta2,0,\
(a31*(2*a2*a2+a32*a32)+a32*(a1*a1-a2*a2))/delta2,2*a1*a2*(a31-a32)/delta2,0},\
{a2*(2*a1*a1-a32*a32)/(a1*delta1),2*(a1*a1*a31+a32*(a2*a2+s))/delta2,0},\
{0,0,1/delta1,0,0,1/delta1,0,0,1/delta1}};
complex N[9][9]={{1,0,0,1/2*((a1*a1*x12)+(a2*a2+a31*a32)*y12),a1*a2*(x12-y12)/2,\
a1*a32*(x12-y12)/2,1/2*((a1*a1*x13)+(a2*a2+a31*a32)*y13),a1*a2*(x13-y13)/2,\
a1*a32*(x13-y13)/2},\
{0,1,0,a1*a2*(x12-y12)/2,1/2*(a2*a2*x12+(a1*a1+a31*a32)*y12),a2*a32*(x12-y12)/2,\
a1*a2*(x13-y13)/2,1/2*(a2*a2*x13+(a1*a1+a31*a32)*y13),a2*a32*(x13-y13)/2},\
{0,0,1,a1*a31*(x12-y12)/2,a2*a31*(x12-y12)/2,(a31*a32*x12+(a1*a1+a2*a2)*y12)/2,\
a1*a31*(x12-y12)/2,a2*a31*(x12-y12)/2,(a31*a32*x12+(a1*a1+a2*a2)*y12)/2},\
{(a1*a1*x21+(a2*a2+a31*a32)*y21)/2,a1*a2*(x21-y21)/2,a1*a32*(y21-x21)/2,1,0,0,\
1/2*((a1*a1*x23)+(a2*a2+a31*a32)*y23),a1*a2*(x23-y23)/2,a1*a32*(x12-y23)/2},\
{a1*a2*(x21-y21)/2,(a2*a2*x21+(a1*a1+a31*a32)*y21)/2,a2*a32*(y21-x21)/2,0,1,0,\
a1*a2*(x23-y23)/2,1/2*(a2*a2*x23+(a1*a1+a31*a32)*y23),a2*a32*(x23-y23)/2},\
{a1*a31*(y21-x21)/2,a2*a31*(y21-x21)/2,(a31*a32*x21+(a1*a1+a2*a2)*y21)/2,0,0,1,\
a1*a31*(x23-y23)/2,a2*a31*(x23-y23)/2,(a31*a32*x23+(a1*a1+a2*a2)*y23)/2},\
{(a1*a1*x31+(a2*a2+a31*a32)*y31)/2,a1*a2*(x31-y31)/2,a1*a32*(y31-x31)/2,\
(a1*a1*x32+(a2*a2+a31*a32)*y32)/2,a1*a2*(x32-y32)/2,a1*a32*(y32-x32)/2,1,0,0},\
{a1*a2*(x31-y31)/2,(a2*a2*x31+(a1*a1+a31*a32)*y31)/2,a2*a32*(y31-x31)/2,a1*a2*(x32-y32)/2,(a2*a2*x32+(a1*a1+a31*a32)*y32)/2 ,a2*a32*(y32-x32)/2,0,1,0},{a1*a31*(y31-x31)/2,a2*a31*(y31-x31)/2,(a31*a32*x31+(a1*a1+a2*a2)*y31)/2,\
a1*a31*(y32-x32)/2,a2*a31*(y32-x32)/2,(a31*a32*x32+(a1*a1+a2*a2)*y32)/2,0,0,1.0}};
complex C[9][9];
for (int t=0;t<9;t++){
for (int r=0;r<9;r++){
for (int q=0;q<9;++q){
C[t][r]+=M[t][q]*N[q][r];} } }
det1[0]=C[0][0];
int k=0;
int u=9;
int j;
int i;
for (k=0;k<u-1;k++)
{ for (i=k;i<u;i++)
{ for (j=k;j<u;j++)
{ A[i][j]=C[i][j]-C[i][k]*C[k][j]/C[k][k];}
det1[k+1]=A[k][k]*det1[k];
for (i=k;i<u;i++)
{ for (j=k;i<u;i++)
{ C[i][j]=A[i][j]; } } }
complex d1=det1[u];
cout<<"determinant="<<cabs(real(d1),imag(d1))<<endl;
return d1;
}
struct element
{float zn;};
ofstream& operator <<(ofstream& out,element el)
{ out<<" "<<el.zn<<endl;
return out;}
double main()
{ double h1,h2,h3,amega=20.0,alpha1=0.0,alpha2=0.0,c1=1.5,c2=5.6,x1=0.0,x2=0.0,\
s1=0.0,h_alpha1=0.01,h_alpha2=0.01,h_amega=0.01;
complex func1,func2;
complex f(double,double,double,double);
double s(double,double,double);
double alpha31(double,double,double);
double cabs(double,double);
char ans;
do{
ofstream file1("data.txt",ios::trunc);
cout<<endl<<"h1=";
cin>>h1;
cout<<endl;
cout<<"h2=";
cin>>h2;
cout<<endl<<"h3=";
cin>>h3;
alpha1=h_alpha1;
alpha2=h_alpha2;
x1=amega/c1;
x2=amega/c2;
do{
func1=det(alpha1,alpha2,alpha31(x1,alpha1,alpha2),\
alpha31(x2,alpha1,alpha2),s(x2,alpha2,alpha2),h1,h2,h3);
ofstream file1("data.txt",ios::app);
if (!file1)
{cerr<<"Error open file!";
return 1;}
file1<<alpha1<<" "<<cabs(real(func1),imag(func1))<<endl;
alpha1+=h_alpha1;
alpha2+=h_alpha2; }
while (alpha1<(2*x2));
file1.close();
cout<<" Press 'N' for QUIT ";
cin>>ans; }
while (ans!='N');
}
|