Методы решения жестких и нежестких краевых задач

Эта система решается методом Гаусса с выделением главного элемента.

В точках, расположенных между узлами, решение находиться при помощи решения задач Коши с начальными условиями в i-ом узле:

.

Применять ортонормирование для краевых задач для жестких обыкновенных дифференциальных уравнений оказывается не надо.

8.3 Шпангоут, выражаемый не дифференциальными, а алгебраическими уравнениями

Рассмотрим случай, когда шпангоут (в точке ) выражается не через дифференциальные уравнения, а через алгебраические уравнения.

Выше мы записывали, что:

Можем представить вектор силовых факторов и перемещений в виде:

,

где - вектор перемещений, - вектор сил и моментов.

Алгебраическое уравнение для шпангоута:

,

где G - матрица жесткости шпангоута, R - вектор перемещений шпангоута, - вектор силовых факторов, которые действуют на шпангоут.

В точке шпангоута имеем:

,

то есть нет разрыва в перемещениях , но есть результирующий вектор силовых факторов , который складывается из сил и моментов слева плюс сил и моментов справа от точки шпангоута.



,

,

,

,

, где ,

что справедливо, если мы не забываем, что в данном случае имеем:

,

то есть вектор перемещений и силовых факторов составляется сначала из перемещений (выше) , а потом из силовых факторов (ниже) .

Здесь необходимо вспомнить, что вектор перемещений выражается через искомый вектор состояния :

,

,

где для удобства было введено переобозначение .

Тогда можем записать:

,

Запишем матричные уравнения для этого случая:

,

,

.

Распишем здесь в уравнении вектор :

,

.

Для обеспечения негромоздкости введем обозначение:

.

Тогда уравнение



примет вид:

.

Для удобства переставим слагаемые в матричных уравнениях, чтобы итоговая система линейных алгебраических уравнений записывалась очевидно:

,

,

.

Таким образом, получаем итоговую систему линейных алгебраических уравнений:

Если к шпангоуту приложено внешнее силовое-моментное воздействие , то

следует переписать в виде , тогда:

.

Тогда матричное уравнение



примет вид:

,

.

Итоговая система линейных алгебраических уравнений примет вид:

8.4 Случай, когда уравнения (оболочки и шпангоута) выражаются не через абстрактные вектора, а через вектора, состоящие из конкретных физических параметров

Рассмотрим случай, когда части оболочечной конструкции и шпангоут выражаются через вектора состояния (типа ), которые (в частном случае) совпадают с векторами физических параметров (типа - перемещения, угол, силы, момент). Тогда матрицы типа будут единичными: . И пусть положительные направления физических параметров одинаковы для всех частей оболочки и шпангоута ().

Тогда будем иметь уравнения:

,

,

,

в виде:

,

,

,

где E - единичная матрица.

Уравнения

,

,

,

,

,

примут вид:

,

,

,

,

.

А уравнения

,



.

примут вид:

,

, где

Итоговая система линейных алгебраических уравнений

где .

Это означает, что уравнение

принимает вид:

, (нет скачка в перемещениях и угле) и

- равновесие шпангоута,



то есть:

(перемещения и угол: нет разрыва)

, где (силы и момент: равновесие).



Глава 9. Анализ и упрощение метода А.А. Абрамова

Метод Абрамова является вариантом метода переноса краевых условий краевой задачи для системы "жестких" обыкновенных дифференциальных уравнений, который также получил распространение на практике. Рассмотрим его.

Введем переобозначения для краевых условий, чтобы записать уравнения в том виде, в котором они записаны А.А. Абрамовым [2, 16]:

где Н(0) - прямоугольная матрица условий на крае x=0, r(0) - вектор внешнего воздействия на край x=0.

В методе А.А. Абрамова задача сводится к решению задачи Коши для прогоночных матрицы Н(х) и вектора r(х). Дифференциальные уравнения переноса краевых условий имеют вид [2,16]:

,

.

Начальные условия для интегрирования этих уравнений имеют вид .

В эти уравнения входит получаемая при их выводе произвольная матрица К(х), и, следовательно, задача не имеет единственного решения. Произвол в выборе матрицы К(х) позволяет наложить на прогоночную матрицу Н(х) дополнительные требования с тем, чтобы численное решение задачи Коши было устойчивым.

В методе А.А. Абрамова элементы прогоночной матрицы Н(х) гарантированы от неограниченного роста, так как на эту матрицу накладывается требование

.

где - транспонированная матрица.

Матрица - квадратная невырожденная.

Из постоянства матрицы следует постоянство всех ее элементов и их ограниченность.

В методе А.А. Абрамова матрица К(х) имеет вид

.

Из этого следует, что уравнения переноса краевых условий в методе А.А. Абрамова имеют вид

,

.

Значения вектора состояния в любой точке определяются встречной прогонкой, то есть путем переноса в эту точку граничных условий как слева, так и справа.

Предлагается путем тождественных преобразований краевых условий матрицу

сделать единичной.

Для этого следует применить к краевым условиям на левом крае предлагавшуюся в усовершенствовании метода С.К. Годунова процедуру построчного ортонормирования матричного уравнения краевых условий.

Тогда строки матрицы Н(х), вычисляемой по дифференциальным уравнениям метода А.А. Абрамова, всегда будут оставаться ортонормированнымн.

Это означает, что при изменении аргумента х векторы-строки матрицы Н(х) поворачиваются в пространстве не меняя своей длины и взаимонаправленности. Следовательно, в методе Абрамова осуществляется естественное, то есть следующее из самих дифференциальных уравнений метода, ортонормирование прогоночной матрицы Н(х).

Таким образом правая часть уравнений А.А. Абрамова, может быть упрощена.

Можем выполнить преобразование - построчное ортонормирование уравнения

для получения эквивалентного краевого условия в виде

,

где строки матрицы W(0) ортонормированы.

В результате уравнения прогонки метода А.А. Абрамова принимают более простой вид в правой части

,

,

так как



.

Первое из уравнений метода А.А. Абрамова является независимым, так как не содержит искомого вектора r(х) и может быть проинтегрировано отдельно от второго. Второе из этих уравнений зависит от результатов интегрирования первого уравнения, так как содержит искомую матрицу Н(х).

Независимое уравнение включает в себя матрицу А(х), характеризующую свойства рассматриваемой деформируемой системы. Таким образом, результаты его интегрирования зависят только от матрицы А(х), а также от вида условий на краю х=0, то есть от значения матрицы Н(0).

Таким образом, независимое уравнение описывает только свойства рассматриваемой системы и не включает в себя информацию о внешнем воздействии на систему. Это означает, что результаты его интегрирования не зависят ни от внешнего воздействия на систему, то есть от вектора F(x), ни от величины воздействия на краю х=0, то есть от значения вектора r(0).

От внешнего воздействия F(x) и величины воздействия r(0) на краю х=0 зависят результаты интегрирования только второго уравнения. Таким образом, второе уравнение отражает все имеющееся внешнее воздействие на рассматриваемую систему, за исключением внешнего воздействия на правом краю, которое учитывается при встречной прогонке с этого правого края.

В случаях многовариантных расчетов, то есть таких расчетов, когда для одной и той же рассматриваемой деформируемой системы необходимо получить решения для многих вариантов внешнего воздействия на эту деформируемую систему, простейшим способом могло бы быть решение большого числа отдельных задач - для каждого варианта внешнего воздействия. Однако этот простейший способ требует слишком больших затрат на вычисления.

Для решения таких многовариантных задач следует организовать вычисления следующим образом. Первое из уравнений А.А. Абрамова интегрируется только однажды и используется для всех вариантов внешнего воздействия. Для каждого варианта внешнего воздействия второе уравнение интегрируется с использованием результатов интегрирования первого уравнения. Первое из уравнений потребуется заново проинтегрировать только в том случае, если измениться вид краевых условий. Что значительно сокращает время решения многовариантных по внешнему воздействию задач.



Глава 10. Метод решения краевых задач для обыкновенных дифференциальных уравнений только с четными производными

10.1 Разрешающие уравнения задач только с четными производными

Многие прикладные задачи теории пластин, оболочек и конструкций из них математически моделируются с помощью линейных обыкновенных дифференциальных уравнений, содержащих только четные производные [5,16,60,98,101].

Приведем примеры некоторых из задач, разрешающие уравнения которых содержат только четные производные.

Уравнение изгиба балки постоянной жесткости имеет вид [16]

,

где q(x) - поверхностная распределенная нагрузка, у(х) - искомый прогиб балки.

Задача о перемещениях мембраны имеет разрешающее уравнение [16,101]

,

где q, N - внешнее поперечное давление и внутреннее нормальное усилие, w(x,y) - искомая функция прогибов.

Плоское напряженное состояние пластины, нагруженной контурными силами, лежащими в ее срединной плоскости, описывается уравнением [16,101]



,

где - функция напряжений Эри.

Изгиб изотропных пластин определяется уравнением [16,101]

,

где q - нормальная к поверхности пластины распределенная внешняя нагрузка, D - цилиндрическая жесткость пластины; w(x,y) - искомая функция прогибов.

Состояние трехслойной изгибаемой панели, состоящей из тонких несущих слоев и среднего слоя-заполнителя из сот или пенопласта, описывается уравнением [5]

где D - изгибная жесткость трехслойной панели; E,v,h - модуль Юнга, коэффициент Пуассона и толщина крайних слоев; В - жесткость крайних слоев при растяжении и сжатии; Н - толщина заполнителя, С - жесткость заполнителя при поперечном сдвиге; F(x,y) - искомая разрешающая функция.

Задача поперечного изгиба ортотропной пластинки сводится к уравнению только с четными производными из работы [5].

В настоящее время для расчетов используются различные варианты уравнений цилиндрических оболочек. Известно [98], что для круговых цилиндрических оболочек разрешающие уравнения разных авторов - Власова, Гольденвейзера, Новожилова, Флюгге, Бейларда, Даревского, Тимошенко - Лява, Доннелла - Власова - Лурье содержат только четные производные.

Уравнения в форме Власова содержат только четные производные для изотропных и ортотропных круговых цилиндрических оболочек, нагружаемых силовым воздействием [98]:

1) Полубезмоментная теория,

2) Моментная техническая теория,

3) Общая моментная теория,

4) Общая моментная теория для ортотропной оболочки.

Только четные производные содержатся и в разрешающих уравнениях, описывающих деформацию круговых цилиндрических оболочек при воздействии на оболочку температурного поля, постоянного по ее толщине, или поля, характеризующегося перепадом температуры по толщине оболочки [98]: дифференциальный дискретный ортогональный интегрирование

1) Полубезмоментиая теория,

2) Моментная техническая теория,

3) Общая моментная теория.

Дифференциальные уравнения с частными производными сводятся, например, методом Фурье [16], методом прямых [101] или другими методами разделения переменных к обыкновенным дифференциальным уравнениям, содержащим только четные производные.

10.2 Основы метода

Пусть задача описывается следующей системой линейных обыкновенных дифференциальных уравнений в матричном виде

,

.

Здесь B=const, то есть рассмотрим систему обыкновенных дифференциальных уравнений с постоянными коэффициентами.

Если в качестве начальных выбираются условия при произвольном значении то общее решение задачи Коши (задачи с начальными, а не краевыми условиями) при условии B=const имеет вид [Гантмахер]

В сокращенных обозначениях:

.

Вычисления выполняются по формулам [Гантмахер]:

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

Может быть сформулирована краевая задача. Краевые условия должны быть заданы на значения векторов - на левом крае и - на правом крае.

Часто элементы искомой вектор-функции у(х) уравнения не являются искомыми физическими параметрами краевой задачи. В общем случае физические параметры задачи формируются линейными комбинациями элементов вектор-функций . функции

Тогда краевые условия, задаваемые на отдельные физические параметры задачи или на линейные зависимости физических параметров на краях, могут быть представлены в виде:

,

.

Выражение для вектор-функции можно получить дифференцированием выражения для .

Пользуясь представлением матрицы и матрицы в виде матричных рядов и выполняя непосредственное их дифференцирование нетрудно показать, что

,

.

Интеграл

относится к интегралам, зависящим от параметра, причем с переменными пределами интегрирования и может быть продифференцированным.

Собрав полученные производные, запишем здесь искомое выражение для вектор-функции

,



Объединим формулы для в одну

.

Можем тогда записать

,

.

Можем записать в совокупности с краевыми условиями по аналогии с методом из главы 3:

,

.

Тогда получаем одну итоговую систему линейных алгебраических уравнений:



Таким образом, записана разрешающая система линейных алгебраических уравнений в смысле и по аналогии с методом из Главы 3, но для дифференциальных уравнений, которые содержат только четные производные.

В случае жестких дифференциальных уравнений, содержащих только четные производные, может применяться 1) либо построчное ортонормирование переносимых краевых условий по аналогии с Главой 6, либо 2) может выстраиваться ленточная матрица сопрягаемых участков интервала интегрирования по аналогии с Главой 7; с той разницей, что использовать следует полученные здесь формулы для случая только четных производных в разрешающих дифференциальных уравнениях краевых задач.

Основные результаты и выводы.

Для краевых задач, сформулированных на основе систем линейных обыкновенных дифференциальных уравнений, получены результаты, позволяющие усовершенствовать известные методы решения этих задач, а также предложить новые методы:

- Предложено усовершенствование метода ортогональной прогонки С.К.Годунова для решения краевых задач с жесткими обыкновенными дифференциальными уравнениями, а именно:

- Предложена формула для начала счета методом прогонки С.К.Годунова

- Предложен второй алгоритм для начала счета методом прогонки С.К.Годунова

- Предложена замена метода численного интегрирования Рунге-Кутты в методе прогонки С.К.Годунова

- Предложены матрично-блочные выводы и реализация алгоритмов начала вычислений для метода С.К.Годунова

- Предложены формулы сопряжения частей интервала интегрирования для метода С.К. Годунова

- Выявлены свойства переноса краевых условий в методе С.К. Годунова

- Предложена модификация метода С.К. Годунова

- Предложен метод «переноса краевых условий» (прямой вариант метода) для решения краевых задач с нежесткими обыкновенными дифференциальными уравнениями

- Предложен метод «дополнительных краевых условий» для решения краевых задач с нежесткими обыкновенными дифференциальными уравнениями

- Предложен метод «дополнительных краевых условий» для решения краевых задач с нежесткими обыкновенными дифференциальными уравнениями

- Предложен метод «половины констант» для решения краевых задач с нежесткими обыкновенными дифференциальными уравнениями

- Предложен метод «переноса краевых условий» (пошаговый вариант метода) для решения краевых задач с жесткими обыкновенными дифференциальными уравнениями, а именно:

- Предложен собственно метод «переноса краевых условий» в произвольную точку интервала интегрирования

- Построен вариант метода для случая «жестких» дифференциальных уравнений

- Предложены формулы для вычисления вектора частного решения неоднородной системы дифференциальных уравнений

- Предложен простейший метод решения краевых задач с жесткими обыкновенными дифференциальными уравнениями без ортонормирования - метод «сопряжения участков интервала интегрирования», которые выражены матричными экспонентами

- Предложен метод расчета оболочек составных и со шпангоутами, а именно:

- Предложен вариант записи метода решения жестких краевых задач без ортонормирования - метода «сопряжения участков, выраженных матричными экспонентами» - через положительные направления формул матричного интегрирования дифференциальных уравнений

- Предложены формулы метода для составных оболочек вращения

- Предложены формулы метода для учета шпангоута, выражаемого не дифференциальными, а алгебраическими уравнениями

- Предложены формулы метода для случая, когда уравнения (оболочки и шпангоута) выражаются не через абстрактные вектора, а через вектора, состоящие из конкретных физических параметров

- Предложены анализ и упрощение метода А.А.Абрамова

- Предложен метод решения краевых задач для обыкновенных дифференциальных уравнений только с четными производными

- Предложены коды программ на С++, размещенные в 3 приложениях

- Вычислительные эксперименты подтвердили результаты, полученные теоретически.



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

1.Абовский Н.П., Андреев Н.П., Деруга А.П., Савченков В.И. Численные методы в теории упругости и теории оболочек. Красно.рек: Изд-во Краснояр. университета, 1986. - 383 с.

2.Абрамов А.А. О переносе граничных условий для систем линейных дифференциальных уравнений (вариант метода прогонки)// Журшш вычислительной математики и математической физики, 1961. - T.I. -N3. -С.542-545.

3.Алфугов Н.А,, Попов Б.Г. Использование операторных матриц для расчета трехслойных цилиндрических оболочек, подкрепленных шпангоутами// Изв. АН СССР, Механика твердого тела, 1977. -N 3. - С. 74-80.

4.Алфугов Н.А., Зиновьев П.А., Попов Б.Г. Расчет многослойных пластин и оболочек из композиционных материалов.

М.Машиностроение, 1984. - 263 с.

5.Амбарцумян С.А. Теория анизотропных пластин. М.:Наука, 1987. - 360 с.

6.Андреев Л.В., Ободан Н.И., Лебедев А.Г. Устойчивость оболочек при неосесимметричной деформации. -М.:Наука, 1988. - 208 с.

7.Бабушка И., Вптасек Э., Прагер М. Численные процессы решения дифференциальных уравнений. -М.:Мир, 1969, - 368 с.

8.Бастатский Б.И. Модификация метода Бубнова *Галеркина в задачах теории пологих оболочек// Строительная механика и расчет сооружений. - М., 1985. -N 5. - С. 4-8.

9.Бахвалов Н.С. Численные методы. Анализ, алгебра, обыковенные дифференциальные уравнения. -М.:Мир, 1984. - 496 с.

10.Белман Р. Введение в теорию матриц. -М.:Наука, 1976. - 352 с.

11.Бенерджи П., Баггерфилд Р. Методы граничных элементов в прикладных науках. -М,:Мир, 1984. - 494 с.

12.Березин И.С., Жидков Н.П. Методы вычислений. М.:Физматгиз, 1962. -Т.1. - 464 с.

13.Березин И.С., Жидков Н.П. Методы вычислений. М.:Физматгиз, 1962. -Т.2. - 635 с.

14.Бидерман В.Л. Применение метода прогонки для численного решения задач строительной механики// Изв. АН СССР МТТ. - М., 1967. -N5. - С. 62-66.

15.Бидерман В.Л. Некоторые вычислительные методы решения задач строительной механики, приводимых к обыкновенным дифференциальным уравнениям. - В кн.:Расчеты на прочность. - М.: Машиностроение, 1976. - С. 8-36.

16.Бидерман В.Л. Механика тонкостенных конструкций. -М.: Машиностроение, 1977. - 488 с.

17.Бидерман В.Л. Прикладная теория механических колебаний. -М.: Высшая школа, 1972. - 416 с.

18.Биргер И.А. Некоторые математические методы решения инженерных задач. - М.: Оборонгиз, 1956. - 152 с.

19.Биргер И.А. Стержни, пластинки и оболочки. М,:Наука, 1992.

-392 с.

20.Блехман И.И., Мышкис А .Д., Пановко Я. Г. Прикладная математпка:нредмет, логика, особенности подходов. -КиевгНаукова думка, 1976. - 269 с.

21.Болотин В.В., Новичков Ю.Н. Механика многослойных конструкций. -М.гМашиностроение, 1980. - 375 с.

22.Бондарь B.C., Даншин В.В., Фролов А.Н. Шаговый метод решения задач нелинейного поведения конструкций// Прикладные проблемы прочности и плястичности. Методы решения задач упругости и пластичности. Всесоюзный межвузовский сб. -Горьковский ун-т, 1986. -С. 26-31.

23.Бронштейн И.Н., Семендяев К.А. Справочник по математике. -М.: Наука, 1986. - 544 с.

24.Буриев Т. Применение разностных схем повышенной точности к решению краевых задач оболочек, плит, стержней и балок// Вопросы вычислительной и прикладной математики. - Ташкент, 1984. -N 75. - С. 51-62.

25.Ван Тассел Д. Стиль, разработка, эффективность, отладка и испытание программ. -М.: Мир, 1985. - 332 с,

26.Вайнберг Д.В., Ждан В.З. Матричные алгоритмы в теории оболочек. -К.:КГУ, 1967. - 165 с.

27.Варвак П.М., Варвак Д.П. Метод сеток в задачах расчета строительных конструкций. - М„: Стройиздат, 1977. - 159 с.

28.Вахитов М.Б. Интегрирующие матрицы - аппарат численного решения дифференциальных уравнений строительной механики// Изв. вузов. Авиационная техника, 1966. -N 3. - С. 50-61.



Приложения

Приложение 1

Программа на С++ расчета цилиндрической оболочки - для метода из главы 6

В качестве проверочных задач использовалась схема консольно закрепленных цилиндрической и сферической оболочек с параметрами R/h=50, 100, 200. Длина цилиндрической оболочки рассматривалась L/R=2, а угловые координаты сферической оболочки рассматривались от /4 до 3/4. На свободном крае рассматривалось нормальное к поверхности оболочек погонное усилие, равномерно распределенное в интервале [-/4, /4]. В качестве среды программирования использовалась система Microsoft Visual Studio 2010 (Visual C++).

Первоначально метод был предложен и обсчитывался в кандидатской диссертации А.Ю.Виноградова в 1993-1995 годах. Тогда оказалось, что без использования ортонормирования в рамках метода успешно решаются задачи осесимметрично нагруженных оболочек вращения. Расчеты тогда выполнялись на компьютере поколения 286. Задачи же неосесимметричного нагружения оболочек вращения можно было решать на компьютерах поколения 286 только с применением процедур построчного ортонормирования - как это и предлагалось в рамках метода. Без процедур ортонормирования в неосесимметричных случаях выдавались только ошибочные графики, представлявшие собой хаотично скачущие большие отрицательные и большие положительные значения, например, изгибающего обезразмеренного момента М1.

Современные компьютеры имеют значительно более совершенное внутреннее устройство и более точные внутренние операции с числами, чем это было в 1993-1995 годах. Поэтому было интересно рассмотреть возможность расчета неосесимметрично нагруженных оболочек, например, цилиндров, на современном аппаратном и программном обеспечении в рамках предложенного метода «переноса краевых условий» совсем без использования процедур построчного ортонормирования.

Оказалось, что неосесимметрично нагруженные цилиндры при некоторых параметрах на современных компьютерах уже можно решать в рамках предложенного метода «переноса краевых условий» совсем без применения операций построчного ортонормирования. Это, например, при параметрах цилиндра L/R=2 и R/h=100.

При параметрах цилиндра L/R=2 и R/h=200 все же оказываются необходимыми процедуры ортонормирования. Но на современных персональных компьютерах уже не наблюдаются сплошные скачки значений от больших отрицательных до больших положительных по всему интервалу между краями цилиндра - как это было на компьютерах поколения 286. В частном случае L/R=2 и R/h=200 наблюдаются лишь незначительные скачки в районе максимума изгибающего обезразмеренного момента М1 на левом крае и небольшой скачек обезразмеренного момента М1 на правом крае.

Приводятся графики изгибающего обезразмеренного момента М1:

- слева приводятся графики, полученные при использовании операций построчного ортонормирования на каждом из 100 шагов, на которые разделялся участок между краями цилиндра,

- справа приводятся графики, полученные совсем без применения операций построчного ортонормирования.

Следует сказать, что в качестве расчетной среды использовалась 32-х битная операционная система Windows XP и среда программирования Microsoft Visual Studio 2010 (Visual C++) использовалась в тех же рамках 32-х битной организации операций с числами. Параметры компьютера такие: ноутбук ASUS M51V (CPU Duo T5800).

Компьютеры будут и дальше развиваться такими же темпами как сейчас и это означает, что в самое ближайшее время для подобных расчетов типа расчета неосесимметрично нагруженных оболочек вращения совсем не потребуется применять ортонормирование в рамках предложенного метода «переноса краевых условий», что существенно упрощает программирование метода и увеличивает скорость расчетов не только по сравнению с другими известными методами, но и по сравнению с собственными характеристиками метода «переноса краевых условий» предыдущих лет.

ПРОГРАММА НА С++ (РАСЧЕТ ЦИЛИНДРА):

//from_A_Yu_Vinogradov.cpp: главный файл проекта.

//Решение краевой задачи - цилиндрической оболочки.

//Интервал интегрирования разбит на 100 участков: левый край - точка 0 и правый край - точка 100

#include "stdafx.h"

#include <iostream>

#include <conio.h>

using namespace std;

//Скалярное произведение векторов - i-й строки матрицы А и j-й строки матрицы С.

double mult(double A[8][8], int i, double C[8][8], int j){

double result=0.0;

for(int k=0;k<8;k++){

result+=A[i][k]*C[j][k];

}

return result;

}

//Вычисление нормы вектора, где вектор это i-я строка матрицы А.

double norma(double A[8][8], int i){

double norma_=0.0;

for(int k=0;k<8;k++){

norma_+=A[i][k]*A[i][k];

}

norma_=sqrt(norma_);

return norma_;

}

//Выполнение ортонормирования. Исходная система A*x=b размерности 8х8 приводиться к системе C*x=d, где строки матрицы С ортонормированы.

void orto_norm_8x8(double A[8][8], double b[8], double C[8][8], double d[8]){

double NORM;

double mult0,mult1,mult2,mult3,mult4,mult5,mult6,mult7;

//Получаем 1-ю строку уравнения C*x=d:

NORM=norma(A,0);

for(int k=0;k<8;k++){

C[0][k]=A[0][k]/NORM;

}

d[0]=b[0]/NORM;

//Получаем 2-ю строку уравнения C*x=d:

mult0=mult(A,1,C,0);

for(int k=0;k<8;k++){

C[1][k]=A[1][k]-mult0*C[0][k];

}

NORM=norma(C,1);

for(int k=0;k<8;k++){

C[1][k]/=NORM;

}

d[1]=(b[1]-mult0*d[0])/NORM;

//Получаем 3-ю строку уравнения C*x=d:

mult0=mult(A,2,C,0); mult1=mult(A,2,C,1);

for(int k=0;k<8;k++){

C[2][k]=A[2][k]-mult0*C[0][k]-mult1*C[1][k];

}

NORM=norma(C,2);

for(int k=0;k<8;k++){

C[2][k]/=NORM;

}

d[2]=(b[2]-mult0*d[0]-mult1*d[1])/NORM;

//Получаем 4-ю строку уравнения C*x=d:

mult0=mult(A,3,C,0); mult1=mult(A,3,C,1); mult2=mult(A,3,C,2);

for(int k=0;k<8;k++){

C[3][k]=A[3][k]-mult0*C[0][k]-mult1*C[1][k]-mult2*C[2][k];

}

NORM=norma(C,3);

for(int k=0;k<8;k++){

C[3][k]/=NORM;

}

d[3]=(b[3]-mult0*d[0]-mult1*d[1]-mult2*d[2])/NORM;

//Получаем 5-ю строку уравнения C*x=d:

mult0=mult(A,4,C,0); mult1=mult(A,4,C,1); mult2=mult(A,4,C,2); mult3=mult(A,4,C,3);

for(int k=0;k<8;k++){

C[4][k]=A[4][k]-mult0*C[0][k]-mult1*C[1][k]-mult2*C[2][k]-

mult3*C[3][k];

}

NORM=norma(C,4);

for(int k=0;k<8;k++){

C[4][k]/=NORM;

}

d[4]=(b[4]-mult0*d[0]-mult1*d[1]-mult2*d[2]-mult3*d[3])/NORM;

//Получаем 6-ю строку уравнения C*x=d:

mult0=mult(A,5,C,0); mult1=mult(A,5,C,1); mult2=mult(A,5,C,2); mult3=mult(A,5,C,3); mult4=mult(A,5,C,4);

for(int k=0;k<8;k++){

C[5][k]=A[5][k]-mult0*C[0][k]-mult1*C[1][k]-mult2*C[2][k]-

mult3*C[3][k]-mult4*C[4][k];

}

NORM=norma(C,5);

for(int k=0;k<8;k++){

C[5][k]/=NORM;

}

d[5]=(b[5]-mult0*d[0]-mult1*d[1]-mult2*d[2]-mult3*d[3]-mult4*d[4])/NORM;

//Получаем 7-ю строку уравнения C*x=d:

mult0=mult(A,6,C,0); mult1=mult(A,6,C,1); mult2=mult(A,6,C,2); mult3=mult(A,6,C,3); mult4=mult(A,6,C,4); mult5=mult(A,6,C,5);

for(int k=0;k<8;k++){

C[6][k]=A[6][k]-mult0*C[0][k]-mult1*C[1][k]-mult2*C[2][k]-

mult3*C[3][k]-mult4*C[4][k]-mult5*C[5][k];

}

NORM=norma(C,6);

for(int k=0;k<8;k++){

C[6][k]/=NORM;

}

d[6]=(b[6]-mult0*d[0]-mult1*d[1]-mult2*d[2]-mult3*d[3]-mult4*d[4]-

mult5*d[5])/NORM;

//Получаем 8-ю строку уравнения C*x=d:

mult0=mult(A,7,C,0); mult1=mult(A,7,C,1); mult2=mult(A,7,C,2); mult3=mult(A,7,C,3); mult4=mult(A,7,C,4); mult5=mult(A,7,C,5);

mult6=mult(A,7,C,6);

for(int k=0;k<8;k++){

C[7][k]=A[7][k]-mult0*C[0][k]-mult1*C[1][k]-mult2*C[2][k]-

mult3*C[3][k]-mult4*C[4][k]-mult5*C[5][k]-mult6*C[6][k];

}

NORM=norma(C,7);

for(int k=0;k<8;k++){

C[7][k]/=NORM;

}

d[7]=(b[7]-mult0*d[0]-mult1*d[1]-mult2*d[2]-mult3*d[3]-mult4*d[4]-

mult5*d[5]-mult6*d[6])/NORM;

}

//Выполнение ортонормирования системы A*x=b с прямоугольной матрицей A коэффициентов размерности 4х8.

void orto_norm_4x8(double A[4][8], double b[4], double C[4][8], double d[4]){

double NORM;

double mult0,mult1,mult2,mult3,mult4,mult5,mult6,mult7;

//Получаем 1-ю строку уравнения C*x=d:

NORM=norma(A,0);

for(int k=0;k<8;k++){

C[0][k]=A[0][k]/NORM;

}

d[0]=b[0]/NORM;

//Получаем 2-ю строку уравнения C*x=d:

mult0=mult(A,1,C,0);

for(int k=0;k<8;k++){

C[1][k]=A[1][k]-mult0*C[0][k];

}

NORM=norma(C,1);

for(int k=0;k<8;k++){

C[1][k]/=NORM;

}

d[1]=(b[1]-mult0*d[0])/NORM;

//Получаем 3-ю строку уравнения C*x=d:

mult0=mult(A,2,C,0); mult1=mult(A,2,C,1);

for(int k=0;k<8;k++){

C[2][k]=A[2][k]-mult0*C[0][k]-mult1*C[1][k];

}

NORM=norma(C,2);

for(int k=0;k<8;k++){

C[2][k]/=NORM;

}

d[2]=(b[2]-mult0*d[0]-mult1*d[1])/NORM;

//Получаем 4-ю строку уравнения C*x=d:

mult0=mult(A,3,C,0); mult1=mult(A,3,C,1); mult2=mult(A,3,C,2);

for(int k=0;k<8;k++){

C[3][k]=A[3][k]-mult0*C[0][k]-mult1*C[1][k]-mult2*C[2][k];

}

NORM=norma(C,3);

for(int k=0;k<8;k++){

C[3][k]/=NORM;

}

d[3]=(b[3]-mult0*d[0]-mult1*d[1]-mult2*d[2])/NORM;

}

//Произведение матрицы A1 размерности 4х8 на матрицу А2 размерности 8х8. Получаем матрицу rezult размерности 4х8:

void mat_4x8_on_mat_8x8(double A1[4][8], double A2[8][8], double rezult[4][8]){

for(int i=0;i<4;i++){

for(int j=0;j<8;j++){

rezult[i][j]=0.0;

for(int k=0;k<8;k++){

rezult[i][j]+=A1[i][k]*A2[k][j];

}

}

}

}

//Умножение матрицы A на вектор b и получаем rezult.

void mat_on_vect(double A[8][8], double b[8], double rezult[8]){

for(int i=0;i<8;i++){

rezult[i]=0.0;

for(int k=0;k<8;k++){

rezult[i]+=A[i][k]*b[k];

}

}

}

//Умножение матрицы A размерности 4х8 на вектор b размерности 8 и получаем rezult размерности 4.

void mat_4x8_on_vect_8(double A[4][8], double b[8], double rezult[4]){

for(int i=0;i<4;i++){

rezult[i]=0.0;

for(int k=0;k<8;k++){

rezult[i]+=A[i][k]*b[k];

}

}

}

//Вычитание из вектора u1 вектора u2 и получение вектора rez=u1-u2. Все вектора размерности 4.

void minus(double u1[4], double u2[4], double rez[4]){

for(int i=0;i<4;i++){

rez[i]=u1[i]-u2[i];

}

}

//Вычисление матричной экспоненты EXP=exp(A*delta_x)

void exponent(double A[8][8], double delta_x, double EXP[8][8]) {

//n - количество членов ряда в экспоненте, m - счетчик членов ряда (m<=n)

int n=100, m;

double E[8][8]={0}, TMP1[8][8], TMP2[8][8];

int i,j,k;

//E - единичная матрица - первый член ряда экспоненты

E[0][0]=1.0; E[1][1]=1.0; E[2][2]=1.0; E[3][3]=1.0;

E[4][4]=1.0; E[5][5]=1.0; E[6][6]=1.0; E[7][7]=1.0;

//первоначальное заполнение вспомогательного массива TMP1 - предыдущего члена ряда для следующего перемножения

//и первоначальное заполнение экспоненты первым членом ряда

for(i=0;i<8;i++) {

for(j=0;j<8;j++) {

TMP1[i][j]=E[i][j];

EXP[i][j]=E[i][j];

}

}

//ряд вычисления экспоненты EXP, начиная со 2-го члена ряда (m=2;m<=n)

for(m=2;m<=n;m++) {

for(i=0;i<8;i++) {

for(j=0;j<8;j++) {

TMP2[i][j]=0;

for(k=0;k<8;k++) {

//TMP2[i][j]+=TMP1[i][k]*A[k][j]*delta_x/(m-1);

TMP2[i][j]+=TMP1[i][k]*A[k][j];

}

TMP2[i][j]*=delta_x;//вынесено за цикл произведения строки на столбец

TMP2[i][j]/=(m-1);//вынесено за цикл произведения строки на столбец

EXP[i][j]+=TMP2[i][j];

}

}

//заполнение вспомогательного массива TMP1 для вычисления следующего члена ряда - TMP2 в следующем шаге цикла по m

if (m<n) {

for(i=0;i<8;i++) {

for(j=0;j<8;j++) {

TMP1[i][j]=TMP2[i][j];

}

}

}

}

}

//Вычисление матрицы MAT_ROW в виде матричного ряда для последующего использования

//при вычислении вектора partial_vector - вектора частного решения неоднородной системы ОДУ на шаге delta_x

void mat_row_for_partial_vector(double A[8][8], double delta_x, double MAT_ROW[8][8]) {

//n - количество членов ряда в MAT_ROW, m - счетчик членов ряда (m<=n)

int n=100, m;

double E[8][8]={0}, TMP1[8][8], TMP2[8][8];

int i,j,k;

//E - единичная матрица - первый член ряда MAT_ROW

E[0][0]=1.0; E[1][1]=1.0; E[2][2]=1.0; E[3][3]=1.0;

E[4][4]=1.0; E[5][5]=1.0; E[6][6]=1.0; E[7][7]=1.0;

//первоначальное заполнение вспомогательного массива TMP1 - предыдущего члена ряда для следующего перемножения

//и первоначальное заполнение MAT_ROW первым членом ряда

for(i=0;i<8;i++) {

for(j=0;j<8;j++) {

TMP1[i][j]=E[i][j];

MAT_ROW[i][j]=E[i][j];

}

}

//ряд вычисления MAT_ROW, начиная со 2-го члена ряда (m=2;m<=n)

for(m=2;m<=n;m++) {

for(i=0;i<8;i++) {

for(j=0;j<8;j++) {

TMP2[i][j]=0;

for(k=0;k<8;k++) {

TMP2[i][j]+=TMP1[i][k]*A[k][j];

}

TMP2[i][j]*=delta_x;

TMP2[i][j]/=m;

MAT_ROW[i][j]+=TMP2[i][j];

}

}

//заполнение вспомогательного массива TMP1 для вычисления следующего члена ряда - TMP2 в следующем шаге цикла по m

if (m<n) {

for(i=0;i<8;i++) {

for(j=0;j<8;j++) {

TMP1[i][j]=TMP2[i][j];

}

}

}

}

}

//Задание вектора внешних воздействий в системе ОДУ - вектора POWER: Y'(x)=A*Y(x)+POWER(x):

void power_vector_for_partial_vector(double x, double POWER[8]){

POWER[0]=0.0;

POWER[1]=0.0;

POWER[2]=0.0;

POWER[3]=0.0;

POWER[4]=0.0;

POWER[5]=0.0;

POWER[6]=0.0;

POWER[7]=0.0;

}

//Вычисление vector - НУЛЕВОГО (частный случай) вектора частного решения

//неоднородной системы дифференциальных уравнений на рассматриваемом участке:

void partial_vector(double vector[8]){

for(int i=0;i<8;i++){

vector[i]=0.0;

}

}

//Вычисление vector - вектора частного решения неоднородной системы дифференциальных уравнений на рассматриваемом участке delta_x:

void partial_vector_real(double expo_[8][8], double mat_row[8][8], double x_, double delta_x, double vector[8]){

double POWER_[8]={0};//Вектор внешней нагрузки на оболочку

double REZ[8]={0};

double REZ_2[8]={0};

power_vector_for_partial_vector(x_, POWER_);//Расчитываем POWER_ при координате x_

mat_on_vect(mat_row, POWER_, REZ);//Умножение матрицы mat_row на вектор POWER_ и получаем вектор REZ

mat_on_vect(expo_, REZ, REZ_2);//Умножение матрицы expo_ на вектор REZ и получаем вектор REZ_2

for(int i=0;i<8;i++){

vector[i]=REZ_2[i]*delta_x;

}

}

//Решение СЛАУ размерности 8 методом Гаусса с выделением главного элемента

int GAUSS(double AA[8][8], double bb[8], double x[8]){

double A[8][8];

double b[8];

for(int i=0;i<8;i++){

b[i]=bb[i];//Работать будем с вектором правых частей b, чтобы исходный вектор bb не изменялся при выходе из подпрограммы

for(int j=0;j<8;j++){

A[i][j]=AA[i][j];//Работать будем с матрицей А, чтобы исходная матрица АА не менялась при выходе из подпрограммы

}

}

int e;//номер строки, где обнаруживается главный (максимальный) коэффициент в столбце jj

double s, t, main;//Вспомогательная величина

for(int jj=0;jj<(8-1);jj++){//Цикл по столбцам jj преобразования матрицы А в верхнетреугольную

e=-1; s=0.0; main=A[jj][jj];

for(int i=jj;i<8;i++){//Находится номер е строки, где лежит главный (максимальный) элемент в столбце jj и делается взаимозамена строк

if ((A[i][jj]*A[i][jj])>s) {//Вместо перемножения (удаляется возможный знак минуса) можно было бы использовать функцию по модулю abs()

e=i; s=A[i][jj]*A[i][jj];

}

}

if (e<0) {

cout<<"Mistake "<<jj<<"\n"; return 0;

}

if (e>jj) {//Если главный элемент не в строке с номером jj. а в строке с номером е

main=A[e][jj];

for(int j=0;j<8;j++){//Взаимная замена двух строк - с номерами e и jj

t=A[jj][j]; A[jj][j]=A[e][j]; A[e][j]=t;

}

t=b[jj]; b[jj]=b[e]; b[e]=t;

}

for(int i=(jj+1);i<8;i++){//Приведение к верхнетреугольной матрице

for(int j=(jj+1);j<8;j++){

A[i][j]=A[i][j]-(1/main)*A[jj][j]*A[i][jj];//Перерасчет коэффициентов строки i>(jj+1)

}

b[i]=b[i]-(1/main)*b[jj]*A[i][jj];

A[i][jj]=0.0;//Обнуляемые элементы столбца под диагональным элементом матрицы А

}

}//Цикл по столбцам jj преобразования матрицы А в верхнетреугольную

x[8-1]=b[8-1]/A[8-1][8-1];//Первоначальное определение последнего элемента искомого решения х (7-го)

for(int i=(8-2);i>=0;i--){//Вычисление елементов решения x[i] от 6-го до 0-го

t=0;

for(int j=1;j<(8-i);j++){

t=t+A[i][i+j]*x[i+j];

}

x[i]=(1/A[i][i])*(b[i]-t);

}

return 0;

}

int main()

{

int nn;//Номер гармоники, начиная с 1-й (без нулевой)

int nn_last=50;//Номер последней гармоники

double Moment[100+1]={0};//Массив физического параметра (момента), что рассчитывается в каждой точке между краями

double step=0.02; //step=(L/R)/100 - величина шага расчета оболочки - шага интервала интегрирования (должна быть больше нуля, т.е. положительная)

double h_div_R;//Величина h/R

h_div_R=1.0/100;

double c2;

c2=h_div_R*h_div_R/12;//Величина h*h/R/R/12

double nju;

nju=0.3;

double gamma;

gamma=3.14159265359/4;//Угол распределения силы по левому краю

//распечатка в файлы:

FILE *fp;

// Open for write

if( (fp = fopen( "C:/test.txt", "w" )) == NULL ) // C4996

printf( "The file 'C:/test.txt' was not opened\n" );

else

printf( "The file 'C:/test.txt' was opened\n" );

for(nn=1;nn<=nn_last;nn++){ //ЦИКЛ ПО ГАРМОНИКАМ, НАЧИНАЯ С 1-ОЙ ГАРМОНИКИ (БЕЗ НУЛЕВОЙ ГАРМОНИКИ)

double x=0.0;//Координата от левого края - нужна для случая неоднородной системы ОДУ для вычисления частного вектора FF

double expo_from_minus_step[8][8]={0};//Матрица для расположения в ней экспоненты на шаге типа (0-x1)

double expo_from_plus_step[8][8]={0};//Матрица для расположения в ней экспоненты на шаге типа (x1-0)

double mat_row_for_minus_expo[8][8]={0};//вспомогательная матрица для расчета частного вектора при движении на шаге типа (0-x1)

double mat_row_for_plus_expo[8][8]={0};//вспомогательная матрица для расчета частного вектора при движении на шаге типа (x1-0)

double MATRIXS[100+1][8][8]={0};//Массив из матриц размерности 8x8 для решения СЛАУ в каждой точке интервала интегрирования

double VECTORS[100+1][8]={0};//Массив векторов правых частей размерности 8 соответствующих СЛАУ

double U[4][8]={0};//Матрица краевых условий левого края размерности 4х8

double u_[4]={0};//Вектор размерности 4 внешнего воздействия для краевых условий левого края

double V[4][8]={0};//Матрица краевых условий правого края размерности 4х8

double v_[4]={0};//Вектор размерности 4 внешнего воздействия для краевых условий правого края

double Y[100+1][8]={0};//Массив векторов-решений соответствующих СЛАУ (в каждой точке интервала между краями): MATRIXS*Y=VECTORS

double A[8][8]={0};//Матрица коэффициентов системы ОДУ

double FF[8]={0};//Вектор частного решения неоднородной ОДУ на участке интервала интегрирования

double Ui[4][8]={0};//Вспомогательная матрица коэффициентов переносимых краевых условий с левого края

double ui_[4]={0};//Правые части переносимых краевых условий с левого края

double ui_2[4]={0};//вспомогательный вектор (промежуточный)

double UiORTO[4][8]={0};//Ортонормированная переносимая матрица с левого края

double ui_ORTO[4]={0};//Соответственно правые части ортонормированного переносимого уравнения с левого края

double Vj[4][8]={0};//Вспомогательная матрица коэффициентов переносимых краевых условий с правого края

double vj_[4]={0};//Правые части переносимых краевых условий с правого края

double vj_2[4]={0};//Вспомогательный вектор (промежуточный)

double VjORTO[4][8]={0};//Ортонормированная переносимая матрица с правого края

double vj_ORTO[4]={0};//Соответственно правые части ортонормированного переносимого уравнения с правого края

double MATRIX_2[8][8]={0};//Вспомогательная матрица

double VECTOR_2[8]={0};//Вспомогательный вектор

double Y_2[8]={0};//Вспомогательный вектор

double nn2,nn3,nn4,nn5,nn6,nn7,nn8;//Возведенный в соответствующие степени номер гармоники nn

nn2=nn*nn; nn3=nn2*nn; nn4=nn2*nn2; nn5=nn4*nn; nn6=nn4*nn2; nn7=nn6*nn; nn8=nn4*nn4;

//Заполнение ненулевых элементов матрицы А коэффициентов системы ОДУ

A[0][1]=1.0;

A[1][0]=(1-nju)/2*nn2; A[1][3]=-(1+nju)/2*nn; A[1][5]=-nju;

A[2][3]=1.0;

A[3][1]=(1+nju)/(1-nju)*nn; A[3][2]=2*nn2/(1-nju); A[3][4]=2*nn/(1-nju);

A[4][5]=1.0;

A[5][6]=1.0;

A[6][7]=1.0;

A[7][1]=-nju/c2; A[7][2]=-nn/c2; A[7][4]=-(nn4+1/c2); A[7][6]=2*nn2;

//Здесь надо первоначально заполнить ненулевыми значениями матрицы и вектора краевых условий U*Y[0]=u_ (слева) и V*Y[100]=v_ (справа) :

U[0][1]=1.0; U[0][2]=nn*nju; U[0][4]=nju; u_[0]=0.0;//Сила T1 на левом крае равна нулю

U[1][0]=-(1-nju)/2*nn; U[1][3]=(1-nju)/2; U[1][5]=(1-nju)*nn*c2; u_[1]=0.0;//Сила S* на левом краю равна нулю

U[2][4]=-nju*nn2; U[2][6]=1.0; u_[2]=0;//Момент M1 на левом краю равен нулю

U[3][5]=(2-nju)*nn2; U[3][7]=-1.0;

u_[3]=-sin(nn*gamma)/(nn*gamma);//Сила Q1* на левом крае распределена на угол -gamma +gamma

V[0][0]=1.0; v_[0]=0.0;//Перемещение u на правом крае равно нулю

V[1][2]=1.0; v_[1]=0.0;//Перемещение v на правом крае равно нулю

V[2][4]=1.0; v_[2]=0.0;//Перемещение w на правом крае равно нулю

V[3][5]=1.0; v_[3]=0.0;//Угол поворота на правом крае равен нулю

//Здесь заканчивается первоначальное заполнение U*Y[0]=u_ и V*Y[100]=v_

orto_norm_4x8(U, u_, UiORTO, ui_ORTO);//Первоначальное ортонормирование краевых условий

orto_norm_4x8(V, v_, VjORTO, vj_ORTO);

//Первоначальное заполнение MATRIXS и VECTORS матричными уравнениями краевых условий соответственно

//UiORTO*Y[0]=ui_ORTO и VjORTO*Y[100]=vj_ORTO:

for(int i=0;i<4;i++){

for(int j=0;j<8;j++){

MATRIXS[0][i][j]=UiORTO[i][j];//Левый край; верхнее матричное уравнение

MATRIXS[100][i+4][j]=VjORTO[i][j];//Правый край (точка номер 101 с индексом 100 - отсчет идет с нуля); нижнее матричное уравнение

}

VECTORS[0][i]=ui_ORTO[i];//Левый край; верхнее матричное уравнение

VECTORS[100][i+4]=vj_ORTO[i];//Правый край (точка номер 101 с индексом 100 - отсчет идет с нуля); нижнее матричное уравнение

}

//Цикл по точкам ii интервала интегрирования заполнения ВЕРХНИХ частей матричных уравнений MATRIXS[ii]*Y[ii]=VECTORS[ii],

//начиная со второй точки - точки с индексом ii=1

exponent(A,(-step),expo_from_minus_step);//Шаг отрицательный (значение шага меньше нуля из-за направления вычисления матричной экспоненты)

x=0.0;//начальное значение координаты - для расчета частного вектора

mat_row_for_partial_vector(A, step, mat_row_for_minus_expo);

for(int ii=1;ii<=100;ii++){

x+=step;//Координата для расчета частного вектора на шаге

mat_4x8_on_mat_8x8(UiORTO,expo_from_minus_step,Ui);//Вычисление матрицы Ui=UiORTO*expo_from_minus_step

//partial_vector(FF);//Вычисление НУЛЕВОГО вектора частного решения системы ОДУ на шаге

partial_vector_real(expo_from_minus_step, mat_row_for_minus_expo, x, (-step),FF);// - для движения слева на право

mat_4x8_on_vect_8(UiORTO,FF,ui_2);//Вычисление вектора ui_2=UiORTO*FF

minus(ui_ORTO, ui_2, ui_);//Вычисление вектора ui_=ui_ORTO-ui_2

orto_norm_4x8(Ui, ui_, UiORTO, ui_ORTO);//Ортонормирование для текущего шага по ii

for(int i=0;i<4;i++){

for(int j=0;j<8;j++){

MATRIXS[ii][i][j]=UiORTO[i][j];

}

VECTORS[ii][i]=ui_ORTO[i];

}

}//Цикл по шагам ii (ВЕРХНЕЕ заполнение)

//Цикл по точкам ii интервала интегрирования заполнения НИЖНИХ частей матричных уравнений MATRIXS[ii]*Y[ii]=VECTORS[ii],

//начиная с предпоследней точки - точки с индексом ii=(100-1) используем ii-- (уменьшение индекса точки)

exponent(A,step,expo_from_plus_step);//Шаг положительный (значение шага больше нуля из-за направления вычисления матричной экспоненты)

x=step*100;//Координата правого края

mat_row_for_partial_vector(A, (-step), mat_row_for_plus_expo);

for(int ii=(100-1);ii>=0;ii--){

x-=step;//Движение справа на лево - для расчета частного вектора

mat_4x8_on_mat_8x8(VjORTO,expo_from_plus_step,Vj);//Вычисление матрицы Vj=VjORTO*expo_from_plus_step

//partial_vector(FF);//Вычисление НУЛЕВОГО вектора частного решения системы ОДУ на шаге

partial_vector_real(expo_from_plus_step, mat_row_for_plus_expo, x, step,FF);// - для движения справа на лево

mat_4x8_on_vect_8(VjORTO,FF,vj_2);//Вычисление вектора vj_2=VjORTO*FF

minus(vj_ORTO, vj_2, vj_);//Вычисление вектора vj_=vj_ORTO-vj_2

orto_norm_4x8(Vj, vj_, VjORTO, vj_ORTO);//Ортонормирование для текущего шага по ii

for(int i=0;i<4;i++){

for(int j=0;j<8;j++){

MATRIXS[ii][i+4][j]=VjORTO[i][j];

}

VECTORS[ii][i+4]=vj_ORTO[i];

}

}//Цикл по шагам ii (НИЖНЕЕ заполнение)

//Решение систем линейных алгебраических уравнений

for(int ii=0;ii<=100;ii++){

for(int i=0;i<8;i++){

for(int j=0;j<8;j++){

MATRIX_2[i][j]=MATRIXS[ii][i][j];//Вспомогательное присвоение для соответствия типов в вызывающей функции GAUSS

}

VECTOR_2[i]=VECTORS[ii][i];//Вспомогательное присвоение для соответствия типов в вызывающей функции GAUSS

}

GAUSS(MATRIX_2,VECTOR_2,Y_2);

for(int i=0;i<8;i++){

Y[ii][i]=Y_2[i];

}

}

//Вычисление момента во всех точках между краями

for(int ii=0;ii<=100;ii++){

Moment[ii]+=Y[ii][4]*(-nju*nn2)+Y[ii][6]*1.0;//Момент M1 в точке [ii]

//U[2][4]=-nju*nn2; U[2][6]=1.0; u_[2]=0;//Момент M1

}

}//ЦИКЛ ПО ГАРМОНИКАМ ЗДЕСЬ ЗАКАНЧИВАЕТСЯ

for(int ii=0;ii<=100;ii++){

fprintf(fp,"%f\n",Moment[ii]);

}

fclose(fp);

printf( "PRESS any key to continue...\n" );

_getch();

return 0;

}

Приложение 2

Программа на С++ расчета сферической оболочки (переменные коэффициенты) - для метода из главы 6.

ПРОГРАММА НА С++ (РАСЧЕТ СФЕРЫ):

//sfera_from_A_Yu_Vinogradov.cpp: главный файл проекта.

//Решение краевой задачи с переменными коэффициентами - сфера.

//Интервал интегрирования разбит на 100 участков: левый край - точка 0 и правый край - точка 100

#include "stdafx.h"

#include <iostream>

#include <conio.h>

#include <math.h> //for tan()

using namespace std;

//Вычисление для гармоники с номером nn для значения переменной (угла) angle_fi - матрицы A_perem 8x8 коэффициентов системы ОДУ

void A_perem_coef(double nju, double c2, int nn, double angle_fi, double A_perem[8][8]){

double nn2,nn3,nn4,nn5,nn6,nn7,nn8;//Возведенный в соответствующие степени номер гармоники nn

nn2=nn*nn; nn3=nn2*nn; nn4=nn2*nn2; nn5=nn4*nn; nn6=nn4*nn2; nn7=nn6*nn; nn8=nn4*nn4;

for(int i=0;i<8;i++){

for(int j=0;j<8;j++){

A_perem[i][j]=0.0;//Первоначальное обнуление матрицы

}

}

//Заполнение ненулевых элементов матрицы А коэффициентов системы ОДУ

A_perem[0][1]=1.0;

A_perem[1][0]=(1-nju)*nn2/2/sin(angle_fi)/sin(angle_fi)+nju+1.0/tan(angle_fi)/tan(angle_fi);

A_perem[1][1]=-1.0/tan(angle_fi);

A_perem[1][2]=(3-nju)/2/sin(angle_fi)/tan(angle_fi);

A_perem[1][3]=-(1+nju)*nn/2/sin(angle_fi);

A_perem[1][5]=-(1+nju);

A_perem[2][3]=1.0;

A_perem[3][0]=(3-nju)*nn/(1-nju)/sin(angle_fi)/tan(angle_fi);

A_perem[3][1]=(1+nju)*nn/(1-nju)/sin(angle_fi);

A_perem[3][2]=2*nn2/(1-nju)/sin(angle_fi)/sin(angle_fi)-1.0+1.0/tan(angle_fi)/tan(angle_fi);

A_perem[3][3]=-1.0/tan(angle_fi);

A_perem[3][4]=(1+nju)*2*nn/(1-nju)/sin(angle_fi);

A_perem[4][5]=1.0;

A_perem[5][6]=1.0;

A_perem[6][7]=1.0;

A_perem[7][0]=-(1+nju)/tan(angle_fi)/c2;

A_perem[7][1]=-(1+nju)/c2;

A_perem[7][2]=-(1+nju)*nn/c2/sin(angle_fi);

A_perem[7][4]=nn2/sin(angle_fi)/sin(angle_fi)*(2+(2-nn2)/sin(angle_fi)/sin(angle_fi)+2.0/tan(angle_fi)/tan(angle_fi))-2*(1+nju)/c2;

A_perem[7][5]=(-2.0-(2*nn2+1)/sin(angle_fi)/sin(angle_fi))/tan(angle_fi);

A_perem[7][6]=-1.0+(2*nn2+1)/sin(angle_fi)/sin(angle_fi);

A_perem[7][7]=-2.0/tan(angle_fi);

}

//Задание вектора внешних воздействий в системе ОДУ - вектора POWER: Y'(x)=A*Y(x)+POWER(x):

void power_vector_for_partial_vector(double x, double POWER[8]){

POWER[0]=0.0;

POWER[1]=0.0;

POWER[2]=0.0;

POWER[3]=0.0;

POWER[4]=0.0;

POWER[5]=0.0;

POWER[6]=0.0;

POWER[7]=0.0;

}

//Скалярное произведение векторов - i-й строки матрицы А и j-й строки матрицы С.

double mult(double A[8][8], int i, double C[8][8], int j){

double result=0.0;

for(int k=0;k<8;k++){

result+=A[i][k]*C[j][k];

}

return result;

}

//Вычисление нормы вектора, где вектор это i-я строка матрицы А.

double norma(double A[8][8], int i){

double norma_=0.0;

for(int k=0;k<8;k++){

norma_+=A[i][k]*A[i][k];

}

norma_=sqrt(norma_);

return norma_;

}

//Выполнение ортонормирования. Исходная система A*x=b размерности 8х8 приводиться к системе C*x=d, где строки матрицы С ортонормированы.

void orto_norm_8x8(double A[8][8], double b[8], double C[8][8], double d[8]){

double NORM;

double mult0,mult1,mult2,mult3,mult4,mult5,mult6,mult7;

//Получаем 1-ю строку уравнения C*x=d:

NORM=norma(A,0);

for(int k=0;k<8;k++){

C[0][k]=A[0][k]/NORM;

}

d[0]=b[0]/NORM;

//Получаем 2-ю строку уравнения C*x=d:

mult0=mult(A,1,C,0);

for(int k=0;k<8;k++){

C[1][k]=A[1][k]-mult0*C[0][k];

}

NORM=norma(C,1);

for(int k=0;k<8;k++){

C[1][k]/=NORM;

}

d[1]=(b[1]-mult0*d[0])/NORM;

//Получаем 3-ю строку уравнения C*x=d:

mult0=mult(A,2,C,0); mult1=mult(A,2,C,1);

for(int k=0;k<8;k++){

C[2][k]=A[2][k]-mult0*C[0][k]-mult1*C[1][k];

}

NORM=norma(C,2);

for(int k=0;k<8;k++){

C[2][k]/=NORM;

}

d[2]=(b[2]-mult0*d[0]-mult1*d[1])/NORM;

//Получаем 4-ю строку уравнения C*x=d:

mult0=mult(A,3,C,0); mult1=mult(A,3,C,1); mult2=mult(A,3,C,2);

for(int k=0;k<8;k++){

C[3][k]=A[3][k]-mult0*C[0][k]-mult1*C[1][k]-mult2*C[2][k];

}

NORM=norma(C,3);

for(int k=0;k<8;k++){

C[3][k]/=NORM;

}

d[3]=(b[3]-mult0*d[0]-mult1*d[1]-mult2*d[2])/NORM;

//Получаем 5-ю строку уравнения C*x=d:

mult0=mult(A,4,C,0); mult1=mult(A,4,C,1); mult2=mult(A,4,C,2); mult3=mult(A,4,C,3);

for(int k=0;k<8;k++){

C[4][k]=A[4][k]-mult0*C[0][k]-mult1*C[1][k]-mult2*C[2][k]-

mult3*C[3][k];

}

NORM=norma(C,4);

for(int k=0;k<8;k++){

C[4][k]/=NORM;

}

d[4]=(b[4]-mult0*d[0]-mult1*d[1]-mult2*d[2]-mult3*d[3])/NORM;

//Получаем 6-ю строку уравнения C*x=d:

mult0=mult(A,5,C,0); mult1=mult(A,5,C,1); mult2=mult(A,5,C,2); mult3=mult(A,5,C,3); mult4=mult(A,5,C,4);

for(int k=0;k<8;k++){

C[5][k]=A[5][k]-mult0*C[0][k]-mult1*C[1][k]-mult2*C[2][k]-

mult3*C[3][k]-mult4*C[4][k];

}

NORM=norma(C,5);

for(int k=0;k<8;k++){

C[5][k]/=NORM;

}

d[5]=(b[5]-mult0*d[0]-mult1*d[1]-mult2*d[2]-mult3*d[3]-mult4*d[4])/NORM;

//Получаем 7-ю строку уравнения C*x=d:

mult0=mult(A,6,C,0); mult1=mult(A,6,C,1); mult2=mult(A,6,C,2); mult3=mult(A,6,C,3); mult4=mult(A,6,C,4); mult5=mult(A,6,C,5);

for(int k=0;k<8;k++){

C[6][k]=A[6][k]-mult0*C[0][k]-mult1*C[1][k]-mult2*C[2][k]-

mult3*C[3][k]-mult4*C[4][k]-mult5*C[5][k];

}

NORM=norma(C,6);

for(int k=0;k<8;k++){

C[6][k]/=NORM;

}

d[6]=(b[6]-mult0*d[0]-mult1*d[1]-mult2*d[2]-mult3*d[3]-mult4*d[4]-

mult5*d[5])/NORM;

//Получаем 8-ю строку уравнения C*x=d:

mult0=mult(A,7,C,0); mult1=mult(A,7,C,1); mult2=mult(A,7,C,2); mult3=mult(A,7,C,3); mult4=mult(A,7,C,4); mult5=mult(A,7,C,5);

mult6=mult(A,7,C,6);

for(int k=0;k<8;k++){

C[7][k]=A[7][k]-mult0*C[0][k]-mult1*C[1][k]-mult2*C[2][k]-

mult3*C[3][k]-mult4*C[4][k]-mult5*C[5][k]-mult6*C[6][k];

}

NORM=norma(C,7);

for(int k=0;k<8;k++){

C[7][k]/=NORM;

}

d[7]=(b[7]-mult0*d[0]-mult1*d[1]-mult2*d[2]-mult3*d[3]-mult4*d[4]-

mult5*d[5]-mult6*d[6])/NORM;

}

//Выполнение ортонормирования системы A*x=b с прямоугольной матрицей A коэффициентов размерности 4х8.

...