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

Алексей Юрьевич Виноградов

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

Предлагаются: Усовершенствование метода ортогональной прогонки С.К. Годунова, 3 метода для нежестких случаев краевых задач, 2 метода для жестких случаев краевых задач, 1 метод расчета оболочек составных и со шпангоутами. По сравнению с монографией «Методы решения жестких и нежестких краевых задач» добавлен материал усовершенствования метода С.К.Годунова, добавлено усовершенствование метода дифференциальной прогонки А.А.Абрамова, добавлен метод для краевых задач для обыкновенных дифференциальных уравнений только с четными производными, добавлено графическое предложение метода численного решения дифференциальных уравнений. Сохранены 3 программы на С++, которые реализуют 2 лучших метода из изложенных.

Математика

Исследования

Вуз: Независимое исследование (нет вуза)

ID: 5c5ecbc67966e10554b5c3df

UUID: b878f3c0-0e96-0137-a2d9-525400005860

Язык: Русский

Опубликовано: почти 6 лет назад

Просмотры: 107

-16.16

Алексей Юрьевич Виноградов

Московский государственный технический университет имени НЭ Баумана (МГТУ им. НЭ Баумана)

Комментировать 0

Рецензировать 0

Скачать - 1,5 МБ

Поделиться работой

МЕЖОТРАСЛЕВОЙ НАУЧНО-ИССЛЕДОВАТЕЛЬСКИЙ ИНСТИТУТ ИНСТИТУЦИОНАЛЬНОГО КОНСАЛТИНГА Виноградов А.Ю. Численные методы решения жестких и нежестких краевых задач Монография Москва 2017

УДК 51(075.8) ББК 22.311я73 В 49 Рекомендовано к публикации ученым советом Межотраслевого научно-исследовательского института институционального консалтинга. Рецензенты: Гамонов Евгений Викторович – доктор физико-математических наук, профессор, старший научный сотрудник SITU IBC Варламов Антон Олегович – кандидат технических наук, доцент, старший научный сотрудник АНОО ДПФО "НИПИ" В 49 Виноградов А.Ю. Численные методы решения жестких и нежестких краевых задач: монография / А.Ю. Виноградов. – Москва: National Research, 2017. 112с. ISBN 978-5-9908927-1-2 Предлагаются: Усовершенствование метода ортогональной прогонки С.К. Годунова, 3 метода для нежестких случаев краевых задач, 2 метода для жестких случаев краевых задач, 1 метод расчета оболочек составных и со шпангоутами. По сравнению с монографией «Методы решения жестких и нежестких краевых задач» добавлен материал усовершенствования метода С.К.Годунова, добавлено усовершенствование метода дифференциальной прогонки А.А.Абрамова, добавлен метод для краевых задач для обыкновенных дифференциальных уравнений только с четными производными, добавлено графическое предложение метода численного решения дифференциальных уравнений. Сохранены 3 программы на С++, которые реализуют 2 лучших метода из изложенных. Публикуется в авторской редакции. ISBN 978-5-9908927-1-2 © А.Ю. Виноградов, 2017

Оглавление Введение ........................................................................................................... 5 Глава 1. Известные формулы теории матриц для обыкновенных дифференциальных уравнений ................................................................... 10 Глава 2. Усовершенствование метода ортогональной прогонки С.К. Годунова для решения краевых задач с жесткими обыкновенными дифференциальными уравнениями ........................................................... 12 2.1. Формула для начала счета методом прогонки С.К. Годунова ......... 12 2.2. Второй алгоритм для начала счета методом прогонки С.К.Годунова ...................................................................................................................... 15 2.3. Замена метода численного интегрирования Рунге-Кутты в методе прогонки С.К.Годунова .............................................................................. 16 2.4. Матрично-блочные выводы и реализация алгоритмов начала вычислений для метода С.К. Годунова .................................................... 16 2.5. Сопряжение частей интервала интегрирования для метода С.К. Годунова ...................................................................................................... 18 2.6. Свойства переноса краевых условий в методе С.К. Годунова ......... 19 2.7. Модификация метода С.К. Годунова ................................................ 20 Глава 3. Метод «переноса краевых условий» (прямой вариант метода) для решения краевых задач с нежесткими обыкновенными дифференциальными уравнениями ...........................................................22 Глава 4. Метод «дополнительных краевых условий» для решения краевых задач с нежесткими обыкновенными дифференциальными уравнениями ..................................................................................................23 Глава 5. Метод «половины констант» для решения краевых задач с нежесткими обыкновенными дифференциальными уравнениями ........ 25 Глава 6. Метод «переноса краевых условий» (пошаговый вариант метода) для решения краевых задач с жесткими обыкновенными дифференциальными уравнениями ...........................................................26 6.1. Метод «переноса краевых условий» в произвольную точку интервала интегрирования .......................................................................26 6.2. Случай «жестких» дифференциальных уравнений ........................ 27 6.3. Формулы для вычисления вектора частного решения неоднородной системы дифференциальных уравнений.................................................29 6.4. Применяемые формулы ортонормирования ...................................32 Глава 7. Простейший метод решения краевых задач с жесткими обыкновенными дифференциальными уравнениями без ортонормирования – метод «сопряжения участков интервала интегрирования», которые выражены матричными экспонентами .......34

Глава 8. Расчет оболочек составных и со шпангоутами простейшим методом «сопряжения участков интервала интегрирования» .................36 8.1. Вариант записи метода решения жестких краевых задач без ортонормирования – метода «сопряжения участков, выраженных матричными экспонентами» - через положительные направления формул матричного интегрирования дифференциальных уравнений 36 8.2. Составные оболочки вращения ......................................................... 37 8.3. Шпангоут, выражаемый не дифференциальными, а алгебраическими уравнениями ................................................................39 8.4. Случай, когда уравнения (оболочки и шпангоута) выражаются не через абстрактные вектора, а через вектора, состоящие из конкретных физических параметров ............................................................................42 Глава 9. Анализ и упрощение метода А.А. Абрамова ................................ 45 Глава 10. Метод решения краевых задач для обыкновенных дифференциальных уравнений только с четными производными ........ 48 10.1. Разрешающие уравнения задач только с четными производными ..................................................................................................................... 48 10.2. Основы метода ...................................................................................49 Приложение 1. Постановка задачи, результаты расчетов и программа на С++ расчета цилиндрической оболочки - для метода из главы 6 ............ 53 Приложение 2. Программа на С++ расчета сферической оболочки (переменные коэффициенты) - для метода из главы 6 ............................. 67 Приложение 3. Постановка задачи, результаты расчетов и программа на С++ расчета цилиндра – для метода из главы 7 ....................................... 80 Приложение 4. Метод главы 7 и программа для этого метода на С++ на английском языке ........................................................................................ 89 Приложение 5. Графическое предложение метода численного решения дифференциальных уравнений ................................................................. 105 Список опубликованных работ .................................................................. 107

Введение Актуальность проблемы: Решение проблемы снижения веса конструкций связано с их усложнением и использованием тонкостенных элементов. Даже простейший вариантный способ конструктивной оптимизации требует параметрических исследований на ЭВМ с использованием численных методов решения краевых задач. Самыми известными среди них являются: ‒ конечно-разностные методы построения приближенных решений дифференциальных уравнений с использованием конечноразностных аппроксимаций производных; ‒ различные модификации метода конечных элементов, метод Бубнова-Галеркина, метод Релея-Ритца, основу которых составляют аппроксимации решений дифференциальных уравнений конечными линейными комбинациями заданных функций: полиномов, тригонометрических функций и т.п.; ‒ методы численного определения интегралов обыкновенных дифференциальных уравнений Рунге-Кутты, Вольтерра, Пикара и т.п. Главным успехом методов конечных разностей и конечных элементов является то, что на их основе построены универсальные алгоритмы и созданы пакеты прикладных программ расчета сложных пространственных силовых конструкций. Построенные вычислительные средства способны выявить поток сил в конструкции и, следовательно, самые напряженные ее элементы. Тем не менее, они требуют неоправданно высоких затрат усилий программиста и мощнейших вычислительных средств, когда ставится задача определения напряжений в местах их концентрации. Наиболее очевидная эффективность методов численного интегрирования обыкновенных дифференциальных уравнений состоит в расчете отдельных частей сложных пространственных конструкций и их отдельных тонкостенных элементов с уточнением напряженнодеформированного состояния в местах его быстрого изменения. Эффективность определяется малыми затратами труда программиста, малыми затратами машинного времени и оперативной памяти ЭВМ. Таким образом, повышение эффективности известных численных методов, построение их модификаций и построение новых методов, является актуальной задачей исследований. Предлагаемая научная новизна состоит в следующем: 5

1. Усовершенствован метод ортогональной прогонки С.К. Годунова, 2. Предложен метод «переноса краевых условий» (прямой вариант метода) для решения краевых задач с нежесткими обыкновенными дифференциальными уравнениями, 3. Предложен метод «дополнительных краевых условий» для решения краевых задач с нежесткими обыкновенными дифференциальными уравнениями, 4. Предложен метод «половины констант» для решения краевых задач с нежесткими обыкновенными дифференциальными уравнениями, 5. Предложен метод «переноса краевых условий» (пошаговый вариант метода) для решения краевых задач с жесткими обыкновенными дифференциальными уравнениями, 6. Предложен простейший метод решения краевых задач с жесткими обыкновенными дифференциальными уравнениями без ортонормирования – метод «сопряжения участков интервала интегрирования», которые выражены матричными экспонентами, 7. Предложен простейший метод расчета оболочек составных и со шпангоутами. 8. Усовершенствован метод дифференциальной прогонки А.А. Абрамова. 9. Предложен метод решения краевых задач для обыкновенных дифференциальных уравнений только с четными производными. 10. Графически предложен метод численного решения дифференциальных уравнений. Некоторые работы, на которых основывается изложение методов, опубликованы совместно с д.ф.-м.н. профессором Ю.И.Виноградовым. Вклад д.ф.-м.н. профессора Ю.И. Виноградова в эти совместные публикации заключался либо 1) в обсуждении результатов проверочных расчетов тех формул и методов, которые предложил А.Ю. Виноградов, либо в том, что 2) в дополнение к методам А.Ю. Виноградова было предложено Ю.И. Виноградовым указание, что матрицы Коши можно вычисять не только в виде матричных экспонент, а дополнительно есть возможность их вычислять в смысле функций Коши-Крылова, используя для этого полученные кем-либо аналитические решения систем дифференциальных уравнений строительной механики пластин и оболочек, либо в том, что 3) Ю.И. Виноградов предложил свою, отличную от формулы А.Ю. Виноградова, формулу вычисления вектора частного 6

решения неоднородной системы обыкновенных дифференциальных уравнений, которая выглядит, однако, гораздо более сложной по сравнению с простой формулой А.Ю. Виноградова. Так же в соавторах отдельных статей указаны еще Ю.А. Гусев и Ю.И. Клюев. Их вклад в материал публикаций состоял в выполнении многовариантных проверочных расчетов в соответствии с формулами, алгоритмами и методами, которые предложил А.Ю Виноградов в своей кандидатской диссертации. Кандидатская физ-мат диссертация А.Ю. Виноградова была защищена в 1996 году. Дополнительно можно сказать, что на основе материла из кандидатской диссертации А.Ю. Виноградова были выполнены еще 2 кандидатских физико-математических диссертации под руководством Ю.И. Виноградова, материал которых состоит в основном в многовариантном применении и в проверке рассчетами того, что было предложено А.Ю. Виноградовым в его кандидатской диссертации. В применении к различным конкретным задачам строительной механики тонкостенных оболочек с выявлением и анализом свойств формул, алгоритмов и методов из кандидатской диссертации А.Ю. Виноградова. Вот данные этих 2 диссертаций: Год: 2008 Петров, Виталий Игоревич «Приведение краевых задач к начальным и исследование концентрации напряжений в тонкостенных конструкциях мультипликативным методом» Ученая cтепень: кандидат физико-математических наук Код cпециальности ВАК: 05.13.18 Специальность: Математическое моделирование, численные методы и комплексы программ. Год: 2003 Гусев, Юрий Алексеевич «Мультипликативные алгоритмы переноса краевых условий в задачах механики деформирования оболочек» Ученая степень: кандидат физико-математических наук Код cпециальности ВАК: 01.02.04 Специальность: Механика деформируемого твердого тела. Кроме того, в соответствии с современными возможностями интернета и при наличии новых диссертаций в открытом доступе было выявлено применение материалов из кандидатской диссертации А.Ю. Виноградова с соответствующими ссылками на соответствующие статьи А.Ю. Виноградова в следующих кандидатских и докторских диссертациях технических и физико-математических наук: Год: 2005 РОССИЙСКАЯ АКАДЕМИЯ НАУК СИБИРСКОЕ ОТДЕЛЕНИЕ Институт вычислительных технологий Юрченко Андрей 7

Васильевич ЧИСЛЕННОЕ РЕШЕНИЕ КРАЕВЫХ ЗАДАЧ УПРУГОГО ДЕФОРМИРОВАНИЯ КОМПОЗИТНЫХ ОБОЛОЧЕК ВРАЩЕНИЯ 05.13.18 – математическое моделирование, численные методы и комплексы программ Диссертация на соискание ученой степени кандидата физико-математических наук Новосибирск. Год: 2003 "Методы и алгоритмы определения напряженнодеформированного состояния тонкостенных подкрепленных конструкций вращения из нелинейно-упругого материала" Автор научной работы: Кочетов, Сергей Николаевич Ученая cтепень: кандидат технических наук Место защиты диссертации: Москва Код cпециальности ВАК: 05.23.17 Специальность: Строительная механика. Год: 1998 «Развитие метода суперэлементов применительно к задачам статики и динамики тонкостенных пространственных систем» Автор научной работы: Чеканин, Александр Васильевич Ученая cтепень: доктор технических наук Место защиты диссертации: Москва Код cпециальности ВАК: 05.23.17 Специальность: Строительная механика. Год: 2005 Автор научной работы: Голушко, Сергей Кузьмич «Прямые и обратные задачи механики упругих композитных пластин и оболочек вращения» Ученая cтепень: доктор физико-математических наук Место защиты диссертации: Новосибирск Код cпециальности ВАК: 01.02.04 Специальность: Механика деформируемого твердого тела. Год: 2003 Автор научной работы: Газизов, Хатиб Шарифзянович «Разработка теории и методов расчета динамики, жесткости и устойчивости составных оболочек вращения» Ученая cтепень: доктор технических наук Место защиты диссертации: Уфа Код cпециальности ВАК: 01.02.04 Специальность: Механика деформируемого твердого тела. Год: 2001 Автор научной работы: Шленов, Алексей Юрьевич «Динамика структурно-неоднородных оболочечных конструкций с учетом упруго-пластических свойств материала» Ученая cтепень: кандидат физико-математических наук Место защиты диссертации: Москва Код cпециальности ВАК: 01.02.04 Специальность: Механика деформируемого твердого тела. Год: 1996 Автор научной работы: Рогов, Анатолий Алексеевич «Динамика трубопровода после разрыва» Ученая cтепень: кандидат физико-математических наук Место защиты диссертации: Москва Код cпециальности ВАК: 01.02.04 Специальность: Механика деформируемого твердого тела. Статьи кандидата физико-математических наук А.Ю. Виноградова опубликованы в таких журналах ВАК как: 8

1. Доклады Академии наук РФ – 2 статьи 2. Механика твердого тела РАН – 2 статьи 3. Журнал вычислительной математики и математической физики РАН – 1 статья 4. Математическое моделирование РАН – 2 статьи 5. Фундаментальные исследования – 1 статья 6. Современные проблемы науки и образования – 1 статья 7. Современные наукоемкие технологии – 1 статья. Всего заимствованные в этой работе формулы решения краевых задач для линейных обыкновенных дифференциальных уравнений были взяты только из 4 источников: 1. Абрамов А.А. О переносе граничных условий для систем линейных дифференциальных уравнений (вариант метода прогонки)// Журнал вычислительной математики и математической физики, 1961. T.I. -N3. -С.542-545. 2. Березин И.С., Жидков Н.П. Методы вычислений. М.:Физматгиз, 1962. -Т.2. - 635 с. 3. Гантмахер Ф.Р. Теория матриц. -М,:Наука, 1988. - 548 с, 4. Годунов С,К. О численном решении краевых задач для систем линейных обыкновенных дифференциальных уравнений// Успехи математических наук, 1961. -Т. 16, вып. 3, (99). – с.171-174. Следует сказать, что объем заимствованных и приведенных в работе формул других авторов составляет 5 страниц, а объем предложенных А.Ю.Виноградовым собственных формул составляет 46 страниц. Плюс написанные А.Ю. Виноградовым программы (коды) на современном и наиболее популярном языке программирования С++ составляют 34 страницы. Часть материалов переведена на английский язык. 9

Глава 1. Известные формулы теории матриц для обыкновенных дифференциальных уравнений Изложение всех методов ведется на примере системы дифференциальных уравнений цилиндрической оболочки ракеты – системы обыкновенных дифференциальных уравнений 8-го порядка (после разделения частных производных методом Фурье). Система линейных обыкновенных дифференциальных уравнений имеет вид: Y ( x)  AY ( x)  F ( x) , где Y (x) – искомая вектор-функция задачи размерности 8х1, Y (x) – производная искомой вектор-функции размерности 8х1, A – квадратная матрица коэффициентов дифференциального уравнения размерности 8х8, F (x ) – вектор-функция внешнего воздействия на систему размерности 8х1. Здесь и далее вектора обозначаем жирным шрифтом вместо черточек над буквами Краевые условия имеют вид: UY (0)  u, VY (1)  v, где Y (0) – значение искомой вектор-функции на левом крае х=0 размерности 8х1, U – прямоугольная горизонтальная матрица коэффициентов краевых условий левого края размерности 4х8, u – вектор внешних воздействий на левый край размерности 4х1, Y (1) – значение искомой вектор-функции на правом крае х=1 размерности 8х1, V – прямоугольная горизонтальная матрица коэффициентов краевых условий правого края размерности 4х8, v – вектор внешних воздействий на правый край размерности 4х1. В случае, когда система дифференциальных уравнений имеет матрицу с постоянными коэффициентами A =const, решение задачи Коши имеет вид [Гантмахер]: x Y ( x)  e A( x x 0)Y ( x0 )  e Ax  e  At F (t )dt , x0 где e A( x  x 0)  E  A( x  x0 )  A ( x  x0 ) / 2! A ( x  x0 ) 3 / 3!..., 2 2 3 где E- это единичная матрица. Матричная экспонента ещё может называться матрицей Коши или матрициантом и может обозначаться в виде: K ( x  x0 )  K ( x  x0 )  e A( x x0) . 10

Тогда решение задачи Коши может быть записано в виде: Y ( x)  K ( x  x0 )Y ( x0 )  Y  ( x  x0 ) , где x Y  ( x  x0 )  e Ax  e  At F (t )dt это вектор частного решения x0 неоднородной системы дифференциальных уравнений. Из теории матриц [Гантмахер] известно свойство перемножаемости матричных экспонент (матриц Коши): K ( xi  x0 )  K ( xi  xi 1 )  K ( xi1  xi 2 )  ... K ( x2  x1 )  K ( x1  x0 ) В случае, когда система дифференциальных уравнений имеет матрицу с переменными коэффициентами A  A(x) , решение задачи Коши предлагается, как это известно, искать при помощи свойства перемножаемости матриц Коши. То есть интервал интегрирования разбивается на малые участки и на малых участках матрицы Коши приближенно вычисляются по формуле для постоянной матрицы в экспоненте. А затем матрицы Коши, вычисленные на малых участках, перемножаются: K ( xi  x0 )  K ( xi  xi 1 )  K ( xi1  xi 2 )  ... K ( x2  x1 )  K ( x1  x0 ) , где матрицы Коши приближенно вычисляются по формуле: K ( xi 1  xi )  e A( xi) xi  exp( A( xi )  xi ) , где xi  xi 1  xi . 11

Глава 2. Усовершенствование метода ортогональной прогонки С.К. Годунова для решения краевых задач с жесткими обыкновенными дифференциальными уравнениями 2.1. Формула для начала счета методом прогонки С.К. Годунова Рассмотрим проблему метода прогонки С.К.Годунова. Предположим, что рассматривается оболочка ракеты. Это тонкостенная труба. Тогда система линейных обыкновенных дифференциальных уравнений будет 8-го порядка, матрица A коэффициентов будет иметь размерность 8х8, искомая вектор-функция Y (x) будет иметь размерность 8х1, а матрицы краевых условий будут прямоугольными горизонтальными размерности 4х8. Тогда в методе прогонки С.К.Годунова для такой задачи решение ищется в следующем виде [Годунов]: Y ( x)  Y1 ( x)c1  Y2 ( x)c2  Y3 ( x)c3  Y4 ( x)c4  Y  ( x) , или можно записать в матричном виде: Y ( x)  Yматрица ( x)c  Y  ( x) , где векторы Y1 ( x),Y2 ( x),Y3 ( x),Y4 ( x) - это линейно независимые векторарешения однородной системы дифференциальных уравнений, а вектор Y  (x) - это вектор частного решения неоднородной системы дифференциальных уравнений. Здесь Yматрица ( x)  Y1 ( x),Y2 ( x),Y3 ( x),Y4 ( x) это матрица размерности 8х4, а c это соответствующий вектор размерности 4х1 из искомых констант c1 , c2 , c3 , c4 . Но вообще-то решение для такой краевой задачи с размерностью 8 (вне рамок метода прогонки С.К.Годунова) может состоять не из 4 линейно независимых векторов Y1 ( x),Y2 ( x),Y3 ( x),Y4 ( x) , а полностью из всех 8 линейно независимых векторов-решений дифференциальных уравнений: однородной системы Y ( x)  Y1 ( x)c1  Y2 ( x)c2  Y3 ( x)c3  Y4 ( x)c4   Y5 ( x)c5  Y6 ( x)c6  Y7 ( x)c7  Y8 ( x)c8  Y  ( x). И как раз трудность и проблема метода прогонки С.К.Годунова и состоит в том, что решение ищется только с половиной возможных векторов и констант и проблема в том, что такое решение с половиной констант должно удовлетворять условиям на левом крае (стартовом для 12

прогонки) при всех возможных значениях констант, чтобы потом найти эти константы из условий на правом крае. То есть в методе прогонки С.К.Годунова есть проблема нахождения таких начальных значений векторов Y1 (0),Y2 (0),Y3 (0),Y4 (0),Y  (0) Y1 ( x),Y2 ( x),Y3 ( x),Y4 ( x),Y  ( x) , чтобы можно было начать прогонку с левого края x =0, то есть чтобы удовлетворялись условия UY (0)  u на левом крае при любых значениях констант c1 , c2 , c3 , c4 . Обычно эта трудность «преодолевается» тем, что дифференциальные уравнения записываются не через функционалы, а через физические параметры и рассматриваются самые простейшие условия на простейшие физические параметры, чтобы начальные значения Y1 (0),Y2 (0),Y3 (0),Y4 (0),Y  (0) можно было угадать. То есть задачи со сложными краевыми условиями так решать нельзя: например, задачи с упругими условиями на краях. Ниже предлагается формула для начала вычислений методом прогонки С.К.Годунова. Выполним построчное ортонормирование матричного уравнения краевых условий на левом крае: UY (0)  u , где матрица U прямоугольная и горизонтальная размерности 4х8. В результате получим эквивалентное уравнение краевых условий на левом крае, но уже с прямоугольной горизонтальной матрицей U орто размерности 4х8, у которой будут 4 ортонормированные строки: U ортоY (0)  uорто , где в результате ортонормирования матрицы преобразован в вектор uорто . U вектор u Как выполнять построчное ортонормирование систем линейных алгебраических уравнений можно посмотреть в [Березин, Жидков]. Дополним прямоугольную горизонтальную матрицу U орто до квадратной невырожденной матрицы W : W U орто M , где матрица M размерности 4х8 должна достраивать матрицу U орто до невырожденной квадратной матрицы W размерности 8х8. В качестве строк матрицы M можно взять те краевые условия, то есть выражения тех физических параметров, которые не входят в параметры левого края или линейно независимы с ними. Это вполне 13

возможно, так как у краевых задач столько линейно независимых физических параметров какова размерность задачи, то есть в данном случае их 8 штук и если 4 заданы на левом крае, то ещё 4 можно взять с правого края. Завершим ортонормирование построенной матрицы W , то есть выполним построчное ортонормирование и получим матрицу Wорто размерности 8х8 с ортонормированными строками: Wорто  U орто M орто . Можем записать, что T Yматрица (0)  (M орто ) транспонированная  M орто . Тогда, подставив в формулу метода прогонки С.К.Годунова, получим: Y (0)  Yматрица (0)c  Y  (0) или T Y (0)  M орто c  Y  (0) . Подставим эту последнюю формулу в краевые условия U ортоY (0)  uорто и получим: T U орто [M орто c  Y  (0)]  uорто . Отсюда получаем, что на левом крае константы c уже не на что не влияют, так как T U орто M орто  0 и остается только найти вектор Y  (0) из выражения: U ортоY  (0)  uорто . Но матрица U орто имеет размерность 4х8 и её надо дополнить до квадратной невырожденной, чтобы найти вектор Y  (0) из решения соответствующей системы линейных алгебраических уравнений: U орто M орто Y  (0)  uорто , 0 где 0 - любой вектор, в том числе вектор из нулей. Отсюда получаем при помощи обратной матрицы: U орто  Y (0)  1 uорто M орто 0 . Тогда формула для начала вычислений методом прогонки С.К.Годунова имеет вид: Y (0)  M T орто c 14 U орто M орто 1 uорто 0 .

Из теории матриц [Гантмахер] известно, что если матрица ортонормирована, то её обратная матрица есть её транспонированная матрица. Тогда последняя формула приобретает вид: T Y (0)  M орто c U орто T uорто M орто T T Y (0)  M орто c  U орто , 0 T M орто uорто 0 , T T T Y (0)  M орто c  U орто uорто  M орто 0, T T Y (0)  M орто c  U орто uорто . T Вектора-столбцы из матрицы M орто и вертикальный вектор-свертка T U орто uорто являются линейно независимыми и удовлетворяют краевому условию UY (0)  u . 2.2. Второй алгоритм для начала счета методом прогонки С.К.Годунова Этот алгоритм требует дополнения матрицы краевых условий U до квадратной невырожденной: U . M Начальные значения Y1 (0),Y2 (0),Y3 (0),Y4 (0),Y  (0) находятся из решения следующих систем линейных алгебраических уравнений: U  u , Y (0)  M 0 1 0 U 0 Yi (0)  , где i  , 0 M i 0 0 1 , 0 0 0 0 , 1 0 0 0 , 0 1 где 0 - вектор из нулей размерности 4х1. Вектора-столбцы Yi (0) и вектор-столбец Y  (0) являются линейно независимыми и участвуя в формировании вектора Y (0) удовлетворяют краевому условию UY (0)  u . 15

2.3. Замена метода численного интегрирования РунгеКутты в методе прогонки С.К.Годунова В методе С.К.Годунова, как показано выше, решение ищется в виде: Y ( x)  Yматрица ( x)c  Y  ( x) . На каждом конкретном участке метода прогонки С.К.Годунова между точками ортогонализации можно вместо метода Рунге-Кутты пользоваться теорией матриц и выполнять расчет через матрицу Коши: Yматрица ( x j )  K ( x j  xi )Yматрица ( xi ) . Так выполнять вычисления быстрее, особенно для дифференциальных уравнений с постоянными коэффициентами, так как в случае постоянных коэффициентов достаточно вычислить один раз матрицу Коши на малом участке и в последующем лишь умножать на эту однажды вычисленную матрицу Коши. И аналогично через теорию матриц можно вычислять и вектор Y  (x) частного решения неоднородной системы дифференциальных уравнений. Или для этого вектора отдельно можно использовать метод Рунге-Кутты, то есть можно комбинировать теорию матриц и метод Рунге-Кутты. 2.4. Матрично-блочные выводы и реализация алгоритмов начала вычислений для метода С.К. Годунова Рассмотрим систему линейных выражающую краевые условия при х=0 алгебраических уравнений, UY (0)  u Пусть имеется построенная квадратная невырожденная матрица G U . U* Аналогично запишем вектор g  u , где введенный вектор u* * u является неизвестным. Запишем систему линейных алгебраических уравнений GY (0)  g или в блочном виде U u Y (0)  * . * U u Отсюда следует, что 16

Y (0)  U U* 1 u  G 1 g  Ng . * u Представим N  T T * . Тогда U Y (0)  * U 1 u u  G 1 g  Ng  T T * *  Tu  T * u* . * u u В тоже время помним, что решение краевой задачи ищется в виде Y ( x)  Yматрица ( x)c  Y  ( x) . Сравнивая Y (0)  T *u*  Tu и Y (0)  Yматрица (0)c  Y  (0) при том, что здесь вектор неизвестных констант это u*  c , то получаем начальные значения векторов для начала интегрирования в методе С.К. Годунова: Yматрица (0)  T  и Y  (0)  Tu . Другой матричный вывод можно изложить в следующем виде. Преобразуем систему U u Y (0)  * * U u путем построчного ортонормирования к эквивалентной системе с ортонормированными строками W w Y (0)  * . * W w Тогда можно записать W Y (0)  * W 1 w W  * * w W T w  WT * w W *T w  W T w  W *T w * . * w Делая сравнение двух выражений: Y (0)  Yматрица (0)c  Y  (0) Y (0)  W *T w*  W T w и из того, что c  w * - вектор неизвестных констант, получаем: Yматрица (0)  W T и Y  (0)  W T w . Заметим, что возможен еще один матрично-блочный вывод формул. Переход от системы U u Y (0)  * * U u к системе 17

W w Y (0)  * * W w можно осуществить еще одним способом, заменив построчное GY (0)  g ортонормирование на следующее ортонормированное разложение матрицы G G T  JL где матрица J имеет ортонормированные столбцы, а матрица L верхнетреугольная. Тогда, учитывая правило транспонирования произведения матриц, можем записать G  LT J T . В результате получим GY (0)  g , LT J T Y (0)  g , J T Y (0)  ( LT ) 1 g . Здесь строки матрицы J T ортонормированные. Сравнивая W w Y (0)  * * W w и J T Y (0)  (LT ) 1 g получаем w W u  JT ,  ( LT ) 1 g  ( LT ) 1 * . * * W w u Таким образом, опять получаем ортонормированные начальные значения искомых вектор-функций решения. 2.5. Сопряжение частей интервала интегрирования для метода С.К. Годунова Для автоматизации вычислительного процесса на всем интервале интегрирования, который составлен для сопряженных оболочек с различными физическими и геометрическими параметрами, деформирование которых описывается различными функциями, необходимо иметь процедуры сопряжения соответствующих функций. В общем случае разрешающие функции различных частей интервала интегрирования задачи не имеют физического смысла, а физические параметры задачи выражаются различным образом через эти функции и их производные. Вместе с этим сопряжение смежных 18

участков должно удовлетворять кинематическим и силовым условиям в точке сопряжения. Решить задачу сопряжения частей интервала интегрирования можно следующим образом. Вектор Р, содержащий физические параметры задачи формируется при помощи матрицы М коэффициентов и искомой вектор-функции Y(x): P  MY (x) где М - квадратная невырожденная матрица. Тогда в точке сопряжения х=х* можем записать выражение M Y- ( x* )  P  M Y ( x* ) , где P - вектор, соответствующий дискретному изменению физических параметров при переходе через точку сопряжения от левого участка к правому; индекс "-" означает "слева от точки сопряжения", а индекс "+" означает "справа от точки сопряжения". В методе прогонки Годунова вектор-функция Y (x) задачи на каждом участке ищется в виде Y ( x)  Yматрица ( x)c  Y  ( x) Предположим, что точка сопряжения не совпадает с точкой ортогонального преобразования. Тогда выражение условий сопряжения смежных участков M Y- ( x* )  P  M Y ( x* ) примет вид M  ( Y- матрица ( x)c -  Y   ( x))  P  M  ( Y матрица ( x)c   Y   ( x)) . Если теперь потребовать c -  c то при прямом ходе метода прогонки продолжить интегрирование при переходе точки сопряжения слева направо можно по следующим выражениям: 1 Y матрица ( x)  M  M Y- матрица ( x) , 1 Y* ( x)  M  (M Y-* ( x)  P ) . 2.6. Свойства переноса краевых условий в методе С.К. Годунова При решении краевой задачи для системы "жестких" линейных обыкновенных дифференциальных уравнений методом С.К. Годунова говорят, что осуществляется дискретная ортогонализация по методу 19

Грамма-Шмидта применительно к вектор-функциям, образующим многообразие решений данной задачи, с целью преодоления тенденции вырождения этих вектор-функций в линейно зависимые. Вместе с тем, при реализации метода Годунова одновременно происходит и перенос краевых условий от начального для интегрирования края к другому. Покажем свойства этого переноса. Ранее было записано Y (0)  Yматрица (0)c  Y  (0) Y (0)  W *T w*  W T w Yматрица (0)  W T и Y  (0)  W T w . Тогда можно сказать, что: ‒ вектор w*, который неизвестен, является вектором констант с, ‒ в тоже время вектор w* имеет физический смысл неизвестного на краю х=0 внешнего воздействия на деформируемую систему, ‒ матрица W* является матрицей краевых условий, неизвестных на краю х=0. Из сформулированных положений следует, что перенос граничных условий в методе С.К. Годунова имеет следующий смысл. Продолжение интегрирования, начиная с вектора Y  (0)  W T w , означает перенос "свертки" W T w матричного уравнения краевых условий при х=0 к правому краю х=1. Продолжение интегрирования, начиная с векторов в матрице Yматрица (0) означает, что матрица краевых условий W*, которые неизвестны на краю х=0, переносится на край х= 1. Интегрирование дифференциальных уравнений ведется с целью переноса на край х=1 вектора с, а значит вектора w*, который выражает условия неизвестные на краю х=0. Перенос матрицы W* и вектора w* означает, что матричное уравнение W *Y  (0)  w* краевых условий, которые неизвестны на краю х=0, переносится на край х=1. 2.7. Модификация метода С.К. Годунова Решение в методе С.К. Годунова ищется, как это записано выше, в виде формулы Y ( x)  Yматрица ( x)c  Y  ( x) . 20

Можем записать эту формулу в двух вариантах – в одном случае формула удовлетворяет краевым условиям левого края (индекс L), а в другом – условиям на правом крае (индекс R): YL ( x)  Yматрица L ( x)c L  Y  L ( x) , YR ( x)  Y матрица R ( x)c R  Y  R ( x) . В произвольной точке имеем YL ( x)  YR ( x) . Тогда получаем Yматрица L ( x)c L  Y  L ( x)  Yматрица R ( x)c R  Y  R ( x) , Yматрица L ( x)c L  Y матрица R ( x)c R  Y  R ( x)  Y  L ( x) , Yматрица L ( x)  Yматрица R ( x) cL  Y  R ( x)  Y  L ( x ) . cR То есть получена система линейных алгебраических уравнений традиционного вида с квадратной матрицей коэффициентов Yматрица L ( x)  Yматрица R ( x) для встречного вычисления векторов констант cL . cR 21

Глава 3. Метод «переноса краевых условий» (прямой вариант метода) для решения краевых задач с нежесткими обыкновенными дифференциальными уравнениями Предлагается выполнять интегрирование по формулам теории матриц [Гантмахер] сразу от некоторой внутренней точки интервала интегрирования к краям: Y (0)  K (0  x)Y ( x)  Y  (0  x) , Y (1)  K (1  x)Y ( x)  Y  (1  x) . Подставим формулу для Y (0) в краевые условия левого края и получим: UY (0)  u , U[K (0  x)Y ( x)  Y  (0  x)]  u , UK (0  x)Y ( x)  u - UY  (0  x) . Аналогично для правых краевых условий получаем: VY (1)  v , V [K (1  x)Y ( x)  Y  (1  x)]  v , VK(1  x)Y ( x)  v  VY  (1  x) . То есть получаем два матричных уравнения краевых условий, перенесенные в рассматриваемую точку x : [UK (0  x)]  Y ( x)  u - UY  (0  x) , [VK(1  x)]  Y ( x)  v  VY  (1  x) . Эти уравнения аналогично объединяются в одну систему линейных алгебраических уравнений с квадратной матрицей коэффициентов для нахождения решения Y (x) в любой рассматриваемой точке x : UK (0  x) u  UY  (0  x)  Y ( x)  . VK(1  x) v  VY  (1  x) 22

Глава 4. Метод «дополнительных краевых условий» для решения краевых задач с нежесткими обыкновенными дифференциальными уравнениями Запишем на левом крае ещё одно уравнение краевых условий: MY (0)  m . В качестве строк матрицы M можно взять те краевые условия, то есть выражения тех физических параметров, которые не входят в параметры краевых условий левого края U или линейно независимы с ними. Это вполне возможно, так как у краевых задач столько независимых физических параметров какова размерность задачи, а в параметры краевых условий входит только половина физических параметров задачи. То есть, например, если рассматривается задача об оболочке ракеты, то на левом крае могут быть заданы 4 перемещения. Тогда для матрицы M можно взять параметры сил и моментов, которых тоже 4, так как полная размерность такой задачи – 8. Вектор m правой части неизвестен и его надо найти и тогда можно считать, что краевая задача решена, то есть сведена к задаче Коши, то есть найден вектор Y (0) из выражения: U u , Y (0)  M m то есть вектор Y (0) находится из решения системы линейных алгебраических уравнений с квадратной невырожденной матрицей коэффициентов, состоящей из блоков U и M . Аналогично запишем на правом крае ещё одно уравнение краевых условий: NY (1)  n , где матрица N записывается из тех же соображений дополнительных линейно независимых параметров на правом крае, а вектор n неизвестен. Для правого края тоже справедлива соответствующая система уравнений: V v . Y (1)  N n Запишем Y (1)  K (1  0)Y (0)  Y  (1  0) и подставим в последнюю систему линейных алгебраических уравнений: V v  [K (1  0)Y (0)  Y  (1  0)]  , N n 23

V v V  K (1  0)Y (0)    Y  (1  0) , N n N V v  V  Y  (1  0)  K (1  0)Y (0)  , N n  N  Y  (1  0) V s .  K (1  0)Y (0)  N t Запишем вектор Y (0) через обратную матрицу: U Y (0)  M 1 u m  и подставим в предыдущую формулу: V U  K (1  0) N M 1  u s  m t Таким образом, мы получили систему уравнений вида: B u s ,  m t где матрица B известна, векторы u и s известны, а векторы m и t неизвестны. Разобьем матрицу B на естественные для нашего случая 4 блока и получим: B11 B21 B12 u s   , B22 m t откуда можем записать, что B11u  B12 m  s, B21u  B22 m  t. Следовательно, искомый вектор m вычисляется по формуле: m  B121 ( s  B11u) А искомый вектор n вычисляется через вектор t : t  B21u  B22 m , n  t  N  Y  (1  0) . 24

Глава 5. Метод «половины констант» для решения краевых задач с нежесткими обыкновенными дифференциальными уравнениями В данном методе используется предложенная С.К. Годуновым идея искать решение в виде только с половиной возможных искомых констант, однако и формула для возможности начала такого расчета и дальнейшее применение матричных экспонент (матриц Коши) предложены А.Ю. Виноградовым. Формула для начала счета с левого края только с половиной возможных констант: T T Y (0)  U орто uорто  M орто c, T Y (0)  U орто uорто T M орто  c . Таким образом, записана в матричном виде формула для начала счета с левого края, когда на левом крае удовлетворены краевые условия. Далее запишем VY (1)  v и Y (1)  K (1  0)Y (0)  Y  (1  0) совместно: V [K (1  0)Y (0)  Y  (1  0)]  v , VK(1  0)Y (0)  v  VY  (1  0) и подставим в эту формулу выражение для Y(0): T VK(1  0) U орто uорто T M орто  v  VY  (1  0) c или T VK(1  0) U орто T M орто uорто c  p. Таким образом, мы получили выражение вида: D uорто c  p, где матрица D имеет размерность 4х8 и может быть естественно представлена в виде двух квадратных блоков размерности 4х4: D1 D2  uорто c  p. Тогда можем записать: D1uорто  D2 c  p . Отсюда получаем, что: c  D21 ( p  D1uорто ) . Таким образом, искомые константы найдены. 25

Глава 6. Метод «переноса краевых условий» (пошаговый вариант метода) для решения краевых задач с жесткими обыкновенными дифференциальными уравнениями 6.1. Метод «переноса краевых условий» в произвольную точку интервала интегрирования Полное решение системы дифференциальных уравнений имеет вид Y ( x)  K ( x  x0 )Y ( x0 )  Y  ( x  x0 ) . Или можно записать: Y (0)  K (0  x1 )Y ( x1 )  Y  (0  x1 ) . Подставляем это выражение для Y (0) в краевые условия левого края и получаем: UY (0)  u , U[K (0  x1 )Y ( x1 )  Y  (0  x1 )]  u , UK (0  x1 )Y ( x1 )  u  UY  (0  x1 ) . Или получаем краевые условия, перенесенные в точку x1 : U1Y ( x1 )  u1 , где U1  UK(0  x1 ) и u1  u  UY  (0  x1 ) . Далее запишем аналогично Y ( x1 )  K ( x1  x2 )Y ( x2 )  Y  ( x1  x2 ) И подставим это выражение для Y ( x1 ) в перенесенные краевые условия точки x1 : U1Y ( x1 )  u1 , U1[K ( x1  x2 )Y ( x2 )  Y  ( x1  x2 )]  u1 , U1 K ( x1  x2 )Y ( x2 )  u1  U1Y  ( x1  x2 ) . Или получаем краевые условия, перенесенные в точку x2 : U 2Y ( x2 )  u2 , где U 2  U1 K ( x1  x2 ) и u2  u1  U1Y  ( x1  x2 ) . И так в точку x  переносим матричное краевое условие с левого края и таким же образом переносим матричное краевое условие с правого края. Покажем шаги переноса краевых условий правого края. Можем записать: Y (1)  K (1  xn1 )Y ( xn1 )  Y  (1  xn1 ) Подставляем это выражение для Y (1) в краевые условия правого края и получаем: 26

VY (1)  v , V [ K (1  xn1 )Y ( xn1 )  Y  (1  xn1 )]  v , VK(1  xn1 )Y ( xn1 )  v  VY  (1  xn1 ) Или получаем краевые условия правого края, перенесенные в точку xn1 : Vn1Y ( xn1 )  v n1 , где Vn1  VK(1  xn1 ) и v n1  v  VY  (1  xn1 ) . Далее запишем аналогично Y ( xn1 )  K ( xn1  xn2 )Y ( xn2 )  Y  ( xn1  xn2 ) И подставим это выражение для Y ( xn1 ) в перенесенные краевые условия точки xn1 : Vn1Y ( xn1 )  v n1 , Vn1[ K ( xn1  xn2 )Y ( xn2 )  Y  ( xn1  xn2 )]  v n1 , Vn1 K ( xn1  xn2 )Y ( xn2 )  v n1  Vn1Y  ( xn1  xn2 ) . Или получаем краевые условия, перенесенные в точку xn2 : Vn2Y ( xn2 )  v n2 , где Vn2  Vn1 K ( xn1  xn2 ) и v n2  v n1  Vn1Y  ( xn1  xn2 ) . И так во внутреннюю точку x  интервала интегрирования переносим матричное краевое условие, как показано, и с левого края и таким же образом переносим матричное краевое условие с правого края и получаем: U Y (x )  u , V Y (x  )  v  . Из этих двух матричных уравнений с прямоугольными горизонтальными матрицами коэффициентов очевидно получаем одну систему линейных алгебраических уравнений с квадратной матрицей коэффициентов: U u   Y (x )  . V v 6.2. Случай «жестких» дифференциальных уравнений В случае «жестких» дифференциальных уравнений предлагается применять построчное ортонормирование матричных краевых условий в процессе их переноса в рассматриваемую точку. Для этого формулы 27

ортонормирования систем линейных алгебраических уравнений можно взять в [Березин, Жидков]. То есть, получив U1Y ( x1 )  u1 применяем к этой группе линейных алгебраических уравнений построчное ортонормирование и получаем эквивалентное матричное краевое условие: U1ортоY ( x1 )  u1орто . И теперь уже в это проортонормированное построчно уравнение подставляем Y ( x1 )  K ( x1  x2 )Y ( x2 )  Y  ( x1  x2 ) . И получаем U1орто [ K ( x1  x2 )Y ( x2 )  Y  ( x1  x2 )]  u1орто , U1орто K ( x1  x2 )Y ( x2 )  u1орто  U1ортоY  ( x1  x2 ) . Или получаем краевые условия, перенесенные в точку x2 : U 2Y ( x2 )  u2 , где U 2  U1орто K ( x1  x2 ) и u2  u1орто  U1ортоY  ( x1  x2 ) . Теперь уже к этой группе линейных алгебраических уравнений применяем построчное ортонормирование и получаем эквивалентное матричное краевое условие: U 2ортоY ( x2 )  u2орто И так далее. И аналогично поступаем с промежуточными матричными краевыми условиями, переносимыми с правого края в рассматриваемую точку. В итоге получаем систему линейных алгебраических уравнений с квадратной матрицей коэффициентов, состоящую из двух независимо друг от друга поэтапно проортонормированных матричных краевых условий, которая решается методом Гаусса с выделением главного элемента для получения решения Y ( x ) в рассматриваемой точке x  :  U орто  Vорто  Y (x )  28  uорто  v орто .

6.3. Формулы для вычисления вектора частного решения неоднородной системы дифференциальных уравнений Вместо формулы для вычисления вектора частного решения неоднородной системы дифференциальных уравнений в виде [Гантмахер]:  Y ( x  x0 )  e x Ax e  At F (t )dt x0 предлагается использовать следующую формулу для каждого отдельного участка интервала интегрирования: xj Y  ( x j  xi )  Y  ( x j  xi )  K ( x j  xi )  K ( xi  t ) F (t )dt . xi Правильность приведенной формулы подтверждается следующим: xj Y ( x j  xi )  exp( A( x j  xi ))  exp( A( xi  t ))F (t )dt , xi  xj  Y ( x j  xi )   exp( A( x j  xi )) exp( A( xi  t ))F (t )dt , xi xj  Y ( x j  xi )   exp( A( x j  xi  xi  t ))F (t )dt , xi xj Y  ( x j  xi )   exp( A( x j  t ))F (t )dt , xi xj Y  ( x j  xi )  exp( Ax j )  exp( At) F (t )dt , xi x Y  ( x  xi )  exp( Ax)  exp( At)F (t )dt , xi что и требовалось подтвердить. Вычисление вектора частного решения системы дифференциальных уравнений производиться при помощи представления матрицы Коши под знаком интеграла в виде ряда и интегрирования этого ряда поэлементно: 29

xj Y ( x j  xi )  Y ( x j  xi )  K ( x j  xi )  K ( xi  t ) F (t )dt  xi xj  K ( x j  xi )  ( E  A( xi  t )  A 2 ( xi  t ) 2 / 2!...)F (t )dt  xi xj xj xj 2  K ( x j  xi )(E  F (t )dt  A  ( xi  t ) F (t )dt  A / 2!  ( xi  t ) 2 F (t )dt  ...). xi xi xi   Эта формула справедлива для случая системы дифференциальных уравнений с постоянной матрицей коэффициентов A =const. Рассмотрим вариант, когда шаги интервала интегрирования выбираются достаточно малыми, что позволяет рассматривать вектор F (t ) на участке ( x j  xi ) приближенно в виде постоянной величины F ( xi )  constant, что позволяет вынести этот вектор из под знаков интегралов: xj xj xj Y  ( x j  xi )  K ( x j  xi )(E  dt  A  ( xi  t )dt  A 2 / 2!  ( xi  t ) 2 dt  ...)F (t ). xi xi xi Известно, что при T=(at+b) имеем  T n dt  В нашем случае имеем  (b - t) n dt  xj Тогда получаем  (x xi i  t ) n dt   1 T n 1  const (при n  -1). a(n  1) 1 (b - t) n 1  const (при n  -1). (-1)(n 1) 1 ( xi  x j ) n1 . n 1 Тогда получаем ряд для вычисления вектора частного решения неоднородной системы дифференциальных уравнений на малом участке ( x j  xi ) : Y  ( x j  xi )  K ( x j  xi )  ( E  A( xi  x j ) / 2! A2 ( xi  x j ) 2 / 3!...) ( x j  xi )  F ( xi ). Для случая дифференциальных уравнений с переменными коэффициентами для каждого участка может использоваться осредненная матрица коэффициентов системы Ai  A( xi ) дифференциальных уравнений. Если рассматриваемый участок интервала интегрирования не мал, то предлагаются следующие итерационные (рекуррентные) формулы. Приведем формулы вычисления вектора частного решения, например, Y  ( x3  x0 ) на рассматриваемом участке ( x3  x0 ) через вектора 30

частного решения Y  ( x1  x0 ) , Y  ( x2  x1 ) , Y  ( x3  x2 ) соответствующих подучастков ( x1  x0 ) , ( x2  x1 ) , ( x3  x2 ) . Имеем Y ( x)  K ( x  x0 )Y ( x0 )  Y  ( x  x0 ) . Также имеем формулу для отдельного подучастка: xj Y ( x j  xi )  Y ( x j  xi )  K ( x j  xi )  K ( xi  t ) F (t )dt . xi   Можем записать: Y ( x1 )  K ( x1  x0 )Y ( x0 )  Y  ( x1  x0 ) , Y ( x2 )  K ( x2  x1 )Y ( x1 )  Y  ( x2  x1 ) . Подставим Y ( x1 ) в Y ( x2 ) и получим: Y ( x2 )  K ( x2  x1 )[K ( x1  x0 )Y ( x0 )  Y  ( x1  x0 )]  Y  ( x2  x1 )   K ( x2  x1 ) K ( x1  x0 )Y ( x0 )  K ( x2  x1 )Y  ( x1  x0 )  Y  ( x2  x1 ) . Сравним полученное выражение с формулой: Y ( x2 )  K ( x2  x0 )Y ( x0 )  Y  ( x2  x0 ) и получим, очевидно, что: K ( x2  x0 )  K ( x2  x1 )K ( x1  x0 ) и для частного вектора получаем формулу: Y  ( x2  x0 )  K ( x2  x1 )Y  ( x1  x0 )  Y  ( x2  x1 ) . То есть вектора подучастков Y  ( x1  x0 ),Y  ( x2  x1 ) не просто складываются друг с другом, а с участием матрицы Коши подучастка. Аналогично запишем Y ( x3 )  K ( x3  x2 )Y ( x2 )  Y  ( x3  x2 ) и подставим сюда формулу для Y ( x2 ) и получим: Y ( x3 )  K ( x3  x2 )[K ( x2  x1 ) K ( x1  x0 )Y ( x0 )  K ( x2  x1 )Y  ( x1  x0 )  Y  ( x2  x1 )]   Y  ( x3  x2 )  K ( x3  x2 ) K ( x2  x1 ) K ( x1  x0 )Y ( x0 )   K ( x3  x2 ) K ( x2  x1 )Y  ( x1  x0 )  K ( x3  x2 )Y  ( x2  x1 )  Y  ( x3  x2 ). Сравнив полученное выражение с формулой: Y ( x3 )  K ( x3  x0 )Y ( x0 )  Y  ( x3  x0 ) очевидно, получаем, что: K ( x3  x0 )  K ( x3  x2 )K ( x2  x1 )K ( x1  x0 ) и вместе с этим получаем формулу для частного вектора: Y  ( x3  x0 )  K ( x3  x2 ) K ( x2  x1 )Y  ( x1  x0 )  K ( x3  x2 )Y  ( x2  x1 )  Y  ( x3  x2 ). То есть именно так и вычисляется частный вектор – вектор частного решения неоднородной системы дифференциальных уравнений, то есть так вычисляется, например, частный вектор на Y  ( x3  x0 ) рассматриваемом участке ( x3  x0 ) через вычисленные частные вектора 31

Y  ( x1  x0 ) , Y  ( x2  x1 ) , Y  ( x3  x2 ) соответствующих подучастков ( x1  x0 ) , ( x2  x1 ) , ( x3  x2 ) . 6.4. Применяемые формулы ортонормирования Взято из [Березин, Жидков]. Пусть дана система линейных алгебраических уравнений порядка n: A x =b . Здесь над векторами поставим черточки вместо их обозначения жирным шрифтом. Будем рассматривать строки матрицы A системы как векторы: ai =( a i1 , a i2 ,…, a in ). Ортонормируем эту систему векторов. Первое уравнение системы A x = b делим на n  a12k . k 1 При этом получим: с11 x1 + с12 x 2 +…+ с1n x n = d 1 , c1 =( c11 , c12 ,…, c1n ), где c1k = a1k n ,  a12k k 1 d1 = b1 n  c12k =1. , n k 1  a12k k 1 Второе уравнение системы заменяется на: с 21 x1 + с 22 x 2 +…+ с 2n x n = d 2 , c2 =( c 21, c 22 ,…, c 2n ), / где c 2k = / c 2k , n  c 2/ 2k d 2= d2 , n  c 2/ 2k k 1 k 1 / / c 2k = a 2k -( a 2 , c1 ) c1k , d 2 = b2 -( a 2 , c1 ) d 1 . Аналогично поступаем дальше. Уравнение с номером i примет вид: с i1 x1 + с i2 x 2 +…+ с in x n = d i , ci =( c i1 , c i2 ,…, c in ), / / где c ik = cik n  cik/ 2 k 1 , di= di n ,  cik/ 2 k 1 32

/ c ik = a ik -( ai , c1 ) c1k -( ai , c2 ) c 2k -…-( ai , ci 1 ) ci 1, k , / d i = bi -( ai , c1 ) d 1 -( ai , c2 ) d 2 -…-( ai , ci 1 ) d i 1. Здесь ( ai , c ) – скалярное произведение векторов. j Процесс будет осуществим, если система линейных алгебраических уравнений линейно независима. В результате мы придем к новой системе C x  d , где матрица C будет с ортонормированными строками, то есть обладает свойством C  C T  E , где E - это единичная матрица. 33

Глава 7. Простейший метод решения краевых задач с жесткими обыкновенными дифференциальными уравнениями без ортонормирования – метод «сопряжения участков интервала интегрирования», которые выражены матричными экспонентами Идея преодоления трудностей вычислений путем разделения интервала интегрирования на сопрягаемые участки принадлежит д.ф.м.н. профессору Ю.И.Виноградову (в том числе на этой идее защищена его докторская физ-мат диссертация), а простейшая реализация этой идеи через формулы теории матриц принадлежит к.ф.-м.н. А.Ю.Виноградову. Разделим интервал интегрирования краевой задачи, например, на 3 участка. Будем иметь точки (узлы), включая края: x0 , x1 , x2 , x3 . Имеем краевые условия в виде: UY ( x0 )  u, VY ( x3 )  v. Можем записать матричные уравнения сопряжения участков: Y ( x0 )  K ( x0  x1 )Y ( x1 )  Y  ( x0  x1 ) , Y ( x1 )  K ( x1  x2 )Y ( x2 )  Y  ( x1  x2 ) , Y ( x2 )  K ( x2  x3 )Y ( x3 )  Y  ( x2  x3 ) . Это мы можем переписать в виде, более удобном для нас далее: EY ( x0 )  K ( x0  x1 )Y ( x1 )  Y  ( x0  x1 ) , EY ( x1 )  K ( x1  x2 )Y ( x2 )  Y  ( x1  x2 ) , EY ( x2 )  K ( x2  x3 )Y ( x3 )  Y  ( x2  x3 ) . где E - единичная матрица. Тогда в объединенном матричном виде получаем систему линейных алгебраических уравнений в следующей форме: U 0 0 0 E  K ( x0  x1 ) 0 0 0 E  K ( x1  x2 ) 0 0 0 E  K ( x 2  x3 ) 0 0 0 V u Y ( x0 ) Y  ( x0  x1 ) Y ( x1 )   Y  ( x1  x2 ) . Y ( x2 ) Y  ( x 2  x3 ) Y ( x3 ) v Эта система решается методом Гаусса с выделением главного элемента. В точках, расположенных между узлами, решение находиться при помощи решения задач Коши с начальными условиями в i-ом узле: 34

Y ( x)  K ( x  xi )Y ( xi )  Y  ( x  xi ) . Применять ортонормирование для краевых задач для жестких обыкновенных дифференциальных уравнений оказывается не надо, так как на каждом участке интервала интегрирования вычисление каждой матричной экспоненты выполняется независимо и от начальной единичной (ортонормированной) матрицы, что делает ненужным применение ортонормирования в отличие от метода Годунова, что значительно упрощает программирование по сравнению с методом Годунова. Вычислять матрицы Коши можно не в виде матричных экспонент, а при помощи методов типа Рунге-Кутты от стартовой единичной матрицы, а вектор частного решения неоднородной системы дифференциальных уравнений вычислять на каждом участке методами типа Рунге-Кутты следует от стартового нулевого вектора. В случае применения методов типа Рунге-Кутты оценки погрешностей хорошо известны, что означает, что вычисления можно выполнять с заранее известной точностью. 35

Глава 8. Расчет оболочек составных и со шпангоутами простейшим методом «сопряжения участков интервала интегрирования» 8.1. Вариант записи метода решения жестких краевых задач без ортонормирования – метода «сопряжения участков, выраженных матричными экспонентами» - через положительные направления формул матричного интегрирования дифференциальных уравнений Разделим интервал интегрирования краевой задачи, например, на 3 участка. Будем иметь точки (узлы), включая края: x0 , x1 , x2 , x3 . Имеем краевые условия в виде: UY ( x0 )  u, VY ( x3 )  v. Можем записать матричные уравнения сопряжения участков: Y ( x1 )  K ( x1  x0 )Y ( x0 )  Y  ( x1  x0 ) , Y ( x2 )  K ( x2  x1 )Y ( x1 )  Y  ( x2  x1 ) , Y ( x3 )  K ( x3  x2 )Y ( x2 )  Y  ( x3  x2 ) . Это мы можем переписать в виде, более удобном для нас далее: EY ( x1 )  K ( x1  x0 )Y ( x0 )  Y  ( x1  x0 ) , EY ( x2 )  K ( x2  x1 )Y ( x1 )  Y  ( x2  x1 ) , EY ( x3 )  K ( x3  x2 )Y ( x2 )  Y  ( x3  x2 ) . где E - единичная матрица. В итоге получаем систему линейных алгебраических уравнений: U 0 0 0 u Y ( x0 )   K ( x1  x0 ) E 0 0 Y ( x1  x0 ) Y ( x1 )  Y  ( x2  x1 ) . 0  K ( x2  x1 ) E 0  Y ( x2 ) 0 0  K ( x3  x 2 ) E Y  ( x3  x 2 ) Y ( x3 ) 0 0 0 V v Эта система решается методом Гаусса с выделением главного элемента. Оказывается, что применять ортонормирование не нужно, так как участки интервала интегрирования выбираются такой длинны, что счет на них является устойчивым. В точках вблизи узлов решение находится путем решения соответствующих задач Коши с началом в i-ом узле: Y ( x)  K ( x  xi )Y ( xi )  Y  ( x  xi ) . 36

8.2. Составные оболочки вращения Рассмотрим сопряжения участков составной оболочки вращения. Пусть имеем 3 участка, где каждый участок может выражаться своими дифференциальными уравнениями и физические параметры могут выражаться по-разному – разными формулами на разных участках: В общем случае (на примере участка 12) физические параметры участка (вектор P12 ( x) ) выражаются через искомые параметры системы обыкновенных дифференциальных уравнений этого участка (через вектор Y12 ( x) ) следующим образом: P12 ( x)  M12Y12 ( x) , где матрица M12 - квадратная невырожденная. При переходе точки сопряжения можем записать в общем виде (но на примере точки сопряжения x1 ): P01 ( x1 )  P0112  L0112 P12 ( x1 ) , где P0112 - дискретное приращение физических параметров (сил, моментов) при переходе с участка «01» на участок «12», а матрица L0112 квадратная невырожденная диагональная и состоит из единиц и минус единиц на главной диагонали для установления правильного соответствия принятых положительных направлений сил, моментов, перемещений и углов при переходе с участка «01» на участок «12», которые могут быть разными (в разных дифференциальных уравнениях разных сопрягаемых участков) – в уравнениях слева от точки сопряжения и в уравнениях справа от точки сопряжения. Два последних уравнения при объединении образуют уравнение: M 01Y01 ( x1 )  P0112  L0112 M12Y12 ( x1 ) . В точке сопряжения x2 аналогично получим уравнение: M12Y12 ( x2 )  P1223  L1223 M 23Y23 ( x2 ) . Если бы оболочка состояла бы из одинаковых участков, то мы могли бы записать в объединенном матричном виде систему линейных алгебраических уравнений в следующей форме: 37

U 0 0 0 u Y ( x0 )   K ( x1  x0 ) E 0 0 Y ( x1  x0 ) Y ( x1 )  Y  ( x2  x1 ) . 0  K ( x2  x1 ) E 0  Y ( x2 ) 0 0  K ( x3  x 2 ) E Y  ( x3  x 2 ) Y ( x3 ) 0 0 0 V v Но в нашем случае оболочка состоит из 3 участков, где средний участок можно считать, например, шпангоутом, выражаемым через свои дифференциальные уравнения. Тогда вместо векторов Y ( x0 ) , Y ( x1 ) , Y ( x2 ) , Y ( x3 ) мы должны рассмотреть вектора: Y01 ( x0 ),Y01 ( x1 ),Y12 ( x1 ),Y12 ( x2 ),Y23 ( x2 ),Y23 ( x3 ) . Тогда матричные уравнения UY ( x0 )  u, VY ( x3 )  v. EY ( x1 )  K ( x1  x0 )Y ( x0 )  Y  ( x1  x0 ) , EY ( x2 )  K ( x2  x1 )Y ( x1 )  Y  ( x2  x1 ) , EY ( x3 )  K ( x3  x2 )Y ( x2 )  Y  ( x3  x2 ) примут вид: UY01 ( x0 )  u, VY23 ( x3 )  v. EY01 ( x1 )  K 01 ( x1  x0 )Y01 ( x0 )  Y01* ( x1  x0 ) , M 01Y01 ( x1 )  P0112  L0112 M12Y12 ( x1 ) , EY12 ( x2 )  K12 ( x2  x1 )Y12 ( x1 )  Y12* ( x2  x1 ) , M12Y12 ( x2 )  P1223  L1223 M 23Y23 ( x2 ) , EY23 ( x3 )  K 23 ( x3  x2 )Y23 ( x2 )  Y23* ( x3  x2 ) . После перестановки слагаемых получаем: UY01 ( x0 )  u, VY23 ( x3 )  v.  K 01 ( x1  x0 )Y01 ( x0 )  EY01 ( x1 )  Y01* ( x1  x0 ) , M 01Y01 ( x1 )  L0112 M12Y12 ( x1 )  P0112 ,  K12 ( x2  x1 )Y12 ( x1 )  EY12 ( x2 )  Y12* ( x2  x1 ) , M12Y12 ( x2 )  L1223 M 23Y23 ( x2 )  P1223 ,  K 23 ( x3  x2 )Y23 ( x2 )  EY23 ( x3 )  Y23* ( x3  x2 ) . 38

В итоге мы можем записать алгебраических уравнений: U 0  K 01 ( x1  x0 ) E 0 M 01 0 0 0 0 0 0 0 0 0 0 0 0  L0112 M 12 0  K12 ( x2  x1 ) E 0 M 12 0 0 0 0 итоговую систему линейных 0 0 0 0 0 u Y01 ( x0 ) * 0 Y01 ( x1  x0 ) Y01 ( x1 ) 0  P0112 Y12 ( x1 ) *  Y12 ( x2  x1 ) 0  Y12 ( x2 )  L1223 M 23 0  P1223 Y23 ( x2 ) *  K 23 ( x3  x2 ) E Y23 ( x3  x2 ) Y23 ( x3 ) 0 V v Эта система решается методом Гаусса с выделением главного элемента. В точках, расположенных между узлами, решение находиться при помощи решения задач Коши с начальными условиями в i-ом узле: Y ( x)  K ( x  xi )Y ( xi )  Y  ( x  xi ) . Применять ортонормирование для краевых задач для жестких обыкновенных дифференциальных уравнений оказывается не надо. 8.3. Шпангоут, выражаемый не дифференциальными, а алгебраическими уравнениями Рассмотрим случай, когда шпангоут (в точке x1 ) выражается не через дифференциальные уравнения, а через алгебраические уравнения. Выше мы записывали, что: P01 ( x1 )  P0112  L0112 P12 ( x1 ) Можем представить вектор P01 ( x1 ) силовых факторов и перемещений в виде: P01 ( x1 )  R01 ( x1 ) , S01 ( x1 ) где R01 ( x1 ) - вектор перемещений, S01 ( x1 ) - вектор сил и моментов. Алгебраическое уравнение для шпангоута: GR  S , где G – матрица жесткости шпангоута, R – вектор перемещений шпангоута, S – вектор силовых факторов, которые действуют на шпангоут. 39

В точке шпангоута имеем: R  0, S  GR , то есть нет разрыва в перемещениях R  0 , но есть результирующий вектор силовых факторов S  GR , который складывается из сил и моментов слева плюс сил и моментов справа от точки шпангоута. P01 ( x1 )  R  L0112 P12 ( x1 ) , S P01 ( x1 )  0  L0112 P12 ( x1 ) , GR R01 ( x1 ) R (x ) 0   L0112 12 1 , S 01 ( x1 ) GR S12 ( x1 ) M 01Y01 ( x1 )  0  L0112 M12Y12 ( x1 ) , GR M 01Y01 ( x1 )  g *  L0112 M12Y12 ( x1 ) , где g *  0 , GR что справедливо, если мы не забываем, что в данном случае имеем: P01 ( x1 )  R01 ( x1 ) , S 01 ( x1 ) то есть вектор перемещений и силовых факторов составляется сначала из перемещений (выше) R01 ( x1 ) , а потом из силовых факторов (ниже) S01 ( x1 ) . Здесь необходимо вспомнить, что вектор перемещений R01 ( x1 ) выражается через искомый вектор состояния Y01 ( x1 ) : g*  0 0  , GR GR01 ( x1 ) R01 ( x1 ) M 11p p P01 ( x1 )   M 01Y01 ( x1 )  M Y01 ( x1 )  S 01 ( x1 ) M 21p M 12p  Y01 ( x1 ) , M 22p где для удобства было введено переобозначение M 01  M p . Тогда можем записать: R01 ( x1 )  M 11p g*  0 GR01 ( x1 )  GM p 11 M 12p  Y01 ( x1 ) , 0 0...0 0   Y01 ( x1 ) p p p  Y01 ( x1 )  p M 12  Y01 ( x1 ) G M 11 M 12 G M 11 M 12p Запишем матричные уравнения для этого случая: 40

UY01 ( x0 )  u, VY12 ( x2 )  v. EY01 ( x1 )  K 01 ( x1  x0 )Y01 ( x0 )  Y01* ( x1  x0 ) , M 01Y01 ( x1 )  g *  L0112 M12Y12 ( x1 ) , EY12 ( x2 )  K12 ( x2  x1 )Y12 ( x1 )  Y12* ( x2  x1 ) . Распишем здесь в уравнении вектор g * : 0 M 01Y01 ( x1 )  ( M 01  GM p 11 0 GM p 11  Y01 ( x1 )  L0112 M 12Y12 ( x1 ) , M 12p )  Y01 ( x1 )  L0112 M 12Y12 ( x1 ) . M 12p Для обеспечения негромоздкости введем обозначение: ( M 01  0 GM p 11 M 12p )  M *. Тогда уравнение M 01Y01 ( x1 )  g *  L0112 M12Y12 ( x1 ) примет вид: M *Y01 ( x1 )  L0112 M 12Y12 ( x1 ) . Для удобства переставим слагаемые в матричных уравнениях, чтобы итоговая система линейных алгебраических уравнений записывалась очевидно: UY01 ( x0 )  u, VY12 ( x2 )  v.  K 01 ( x1  x0 )Y01 ( x0 )  EY01 ( x1 )  Y01* ( x1  x0 ) , M *Y01 ( x1 )  L0112 M12Y12 ( x1 )  0 ,  K12 ( x2  x1 )Y12 ( x1 )  EY12 ( x2 )  Y12* ( x2  x1 ) . Таким образом, получаем алгебраических уравнений: U 0  K 01 ( x1  x0 ) E 0 M* 0 0 0 0 Если к шпангоуту воздействие g p , то g*  итоговую систему линейных 0 0 0 u Y01 ( x0 ) * 0 Y01 ( x1  x0 ) Y01 ( x1 ) .   L0112 M 12 0  0 Y12 ( x1 ) *  K12 ( x2  x1 ) E Y12 ( x2  x1 ) Y12 ( x2 ) 0 V v приложено внешнее силовое-моментное 0 0 0  следует переписать в виде g *  , тогда: p GR01 ( x1 )  g p GR  g GR 41

g*  GM p 11 0 0...0 0 . p  p p  Y01 ( x1 )  M  Y01 ( x1 )  g G M 11 M 12 gp p 12 Тогда матричное уравнение M 01Y01 ( x1 )  0 GM p 11 M 12p  Y01 ( x1 )  L0112 M 12Y12 ( x1 ) примет вид: M 01Y01 ( x1 )  0 GM p 11 M  Y01 ( x1 )  p 12 0  L0112 M 12Y12 ( x1 ) , gp M *Y01 ( x1 )  L0112 M 12Y12 ( x1 )   0 . gp Итоговая система линейных алгебраических уравнений примет вид: U 0  K 01 ( x1  x0 ) E 0 M* 0 0 0 0 u 0 * Y01 ( x0 ) Y01 ( x1  x0 ) 0 0 Y (x )  p .  L0112 M 12 0  01 1  g Y12 ( x1 )  K12 ( x2  x1 ) E * Y12 ( x2 ) Y12 ( x2  x1 ) 0 V v 0 0 8.4. Случай, когда уравнения (оболочки и шпангоута) выражаются не через абстрактные вектора, а через вектора, состоящие из конкретных физических параметров Рассмотрим случай, когда части оболочечной конструкции и шпангоут выражаются через вектора состояния (типа Y12 ( x) ), которые (в частном случае) совпадают с векторами физических параметров (типа P12 ( x) - перемещения, угол, силы, момент). Тогда матрицы типа M12 будут единичными: M12  E . И пусть положительные направления физических параметров одинаковы для всех частей оболочки и шпангоута ( L0112  E ). Тогда будем иметь уравнения: P12 ( x)  M12Y12 ( x) , P01 ( x1 )  P0112  L0112 P12 ( x1 ) , M 01Y01 ( x1 )  P0112  L0112 M12Y12 ( x1 ) , в виде: P12 ( x)  EY12 ( x) , P01 ( x1 )  P0112  EP12 ( x1 ) , EY01 ( x1 )  P0112  EY12 ( x1 ) , 42

где E – единичная матрица. Уравнения M 01Y01 ( x1 )  g *  L0112 M12Y12 ( x1 ) , 0  L0112 P12 ( x1 ) , GR01 ( x1 )  g p P01 ( x1 )  R01 ( x1 ) M 11p p P01 ( x1 )   M 01Y01 ( x1 )  M Y01 ( x1 )  S 01 ( x1 ) M 21p R01 ( x1 )  M 11p g*  GM p 11 M 12p  Y01 ( x1 ) , M 22p M 12p  Y01 ( x1 ) , 0 0 0 0 ,  p  p p p  Y01 ( x1 )  M 12  Y01 ( x1 ) G M 11 M 12 g gp примут вид: EY01 ( x1 )  g *  EY12 ( x1 ) , EY01 ( x1 )  P01 ( x1 )  0  EY12 ( x1 ) , GR01 ( x1 )  g p R01 ( x1 ) E 0  EY01 ( x1 )  M pY01 ( x1 )   Y01 ( x1 ) , S 01 ( x1 ) 0 E R01 ( x1 )  E 0  Y01 ( x1 ) , g*  0 GR01 ( x1 )  0 0 0 0 0   p   Y01 ( x1 )  p . p G E 0  Y01 ( x1 ) g GE 0 g g А уравнения M 01Y01 ( x1 )  0 GM p 11 M  Y01 ( x1 )  p 12 0  L0112 M 12Y12 ( x1 ) , gp M *Y01 ( x1 )  L0112 M 12Y12 ( x1 )   0 . gp примут вид: EY01 ( x1 )  0 0  Y01 ( x1 )  p  EY12 ( x1 ) , GE 0 g M *Y01 ( x1 )  EY12 ( x1 )   0 0 )  M* , где ( E  p GE 0 g Итоговая система линейных алгебраических уравнений U 0  K 01 ( x1  x0 ) E 0 M* 0 0 0 0 u 0 * Y01 ( x0 ) Y01 ( x1  x0 ) 0 0 Y (x )  p  L0112 M 12 0  01 1  g Y12 ( x1 )  K12 ( x2  x1 ) E * Y12 ( x2 ) Y12 ( x2  x1 ) 0 V v 0 0 43

примет вид: U 0  K 01 ( x1  x0 ) E 0 M* 0 0 0 0 где M  (8Ex8 0 4 x8 * G E 4x4 4x4 0 u 0 0 Y01 ( x0 ) Y01* ( x1  x0 ) 0 0 0 Y (x )  p , E 0  01 1  g Y12 ( x1 )  K12 ( x2  x1 ) E * Y12 ( x2 ) Y12 ( x2  x1 ) 0 V v )  (E 8 x8 4x4 0 4 x8 G 4x4 0 ) 4x4 E G 0 . E 4x4 4x4 4x4 4x4 Это означает, что уравнение M *Y01 ( x1 )  EY12 ( x1 )   0 gp принимает вид: E G 4x4 4x4 E G 4x4 4x4 0 0 Y01 ( x1 )  EY12 ( x1 )   p E g 4x4 4x4 0 R01 ( x1 ) E  4x4 E S 01 ( x1 ) 0 4x4 4x4 4x4 0 R12 ( x1 ) 0  p E S12 ( x1 ) g 4x4 4x4 R01 ( x1 )  R12 ( x1 )  0 , (нет скачка в перемещениях и угле) и GR01 ( x1 )  S01 ( x1 )  S12 ( x1 )   g p - равновесие шпангоута, то есть: R01 ( x1 )  R12 ( x1 ) (перемещения и угол: нет разрыва) S 01 ( x1 )  S  g p  S12 ( x1 ) , где S  GR01 ( x1 ) (силы и момент: равновесие). 44

Глава 9. Анализ и упрощение метода А.А. Абрамова Метод Абрамова является вариантом метода переноса краевых условий краевой задачи для системы "жестких" обыкновенных дифференциальных уравнений, который также получил распространение на практике. Рассмотрим его. Введем переобозначения для краевых условий, чтобы записать уравнения в том виде, в котором они записаны А.А. Абрамовым: H (0)Y (0)  r (0) где Н(0) - прямоугольная матрица условий на крае x=0, r(0) – вектор внешнего воздействия на край x=0. В методе А.А. Абрамова задача сводится к решению задачи Коши для прогоночных матрицы Н(х) и вектора r(х). Дифференциальные уравнения переноса краевых условий имеют вид: H / ( x)  H ( x) A( x)  K ( x)H ( x) , r / ( x)  H ( x)F ( x)  K ( x)r ( x) . Начальные условия для интегрирования этих уравнений имеют вид H (0), r (0) . В эти уравнения входит получаемая при их выводе произвольная матрица К(х), и, следовательно, задача не имеет единственного решения. Произвол в выборе матрицы К(х) позволяет наложить на прогоночную матрицу Н(х) дополнительные требования с тем, чтобы численное решение задачи Коши было устойчивым. В методе А.А. Абрамова элементы прогоночной матрицы Н(х) гарантированы от неограниченного роста, так как на эту матрицу накладывается требование HH T  H (0)H (0)T  const . где H T - транспонированная матрица. Матрица HH T - квадратная невырожденная. Из постоянства матрицы HH T следует постоянство всех ее элементов и их ограниченность. В методе А.А. Абрамова матрица К(х) имеет вид K ( x)  HAHT (HH T ) 1 . Из этого следует, что уравнения переноса краевых условий в методе А.А.Абрамова имеют вид H /  HAH TH  HA , r /  HAHTr  HF . 45

Значения вектора состояния в любой точке определяются встречной прогонкой, то есть путем переноса в эту точку граничных условий как слева, так и справа. Предлагается путем тождественных преобразований краевых условий матрицу   HH T  H (0)H (0)T  const сделать единичной. Для этого следует применить к краевым условиям на левом крае предлагавшуюся в усовершенствовании метода С.К. Годунова процедуру построчного ортонормирования матричного уравнения краевых условий. Тогда строки матрицы Н(х), вычисляемой по дифференциальным уравнениям метода А.А. Абрамова, всегда будут оставаться ортонормированными. Это означает, что при изменении аргумента х векторы-строки матрицы Н(х) поворачиваются в пространстве не меняя своей длины и взаимонаправленности. Следовательно, в методе Абрамова осуществляется естественное, то есть следующее из самих дифференциальных уравнений метода, ортонормирование прогоночной матрицы Н(х). Таким образом, правая часть уравнений А.А. Абрамова, может быть упрощена. Можем выполнить преобразование – построчное ортонормирование уравнения H (0)Y (0)  r (0) для получения эквивалентного краевого условия в виде W (0)Y (0)  w (0) , где строки матрицы W(0) ортонормированы. В результате уравнения прогонки метода А.А. Абрамова принимают более простой вид в правой части H /  HA(H T H  E) , r / ( x)  HAHT r  HF , так как   HH T  H (0)H (0)T  const  E  единичнаяматрица . Первое из уравнений метода А.А. Абрамова является независимым, так как не содержит искомого вектора r(х) и может быть проинтегрировано отдельно от второго. Второе из этих уравнений зависит от результатов интегрирования первого уравнения, так как содержит искомую матрицу Н(х). 46

Независимое уравнение включает в себя матрицу А(х), характеризующую свойства рассматриваемой деформируемой системы. Таким образом, результаты его интегрирования зависят только от матрицы А(х), а также от вида условий на краю х=0, то есть от значения матрицы Н(0). Таким образом, независимое уравнение описывает только свойства рассматриваемой системы и не включает в себя информацию о внешнем воздействии на систему. Это означает, что результаты его интегрирования не зависят ни от внешнего воздействия на систему, то есть от вектора F(x), ни от величины воздействия на краю х=0, то есть от значения вектора r(0). От внешнего воздействия F(x) и величины воздействия r(0) на краю х=0 зависят результаты интегрирования только второго уравнения. Таким образом, второе уравнение отражает все имеющееся внешнее воздействие на рассматриваемую систему, за исключением внешнего воздействия на правом краю, которое учитывается при встречной прогонке с этого правого края. В случаях многовариантных расчетов, то есть таких расчетов, когда для одной и той же рассматриваемой деформируемой системы необходимо получить решения для многих вариантов внешнего воздействия на эту деформируемую систему, простейшим способом могло бы быть решение большого числа отдельных задач - для каждого варианта внешнего воздействия. Однако этот простейший способ требует слишком больших затрат на вычисления. Для решения таких многовариантных задач следует организовать вычисления следующим образом. Первое из уравнений А.А. Абрамова интегрируется только однажды и используется для всех вариантов внешнего воздействия. Для каждого варианта внешнего воздействия второе уравнение интегрируется с использованием результатов интегрирования первого уравнения. Первое из уравнений потребуется заново проинтегрировать только в том случае, если измениться вид краевых условий. Что значительно сокращает время решения многовариантных по внешнему воздействию задач. 47

Глава 10. Метод решения краевых задач для обыкновенных дифференциальных уравнений только с четными производными 10.1. Разрешающие уравнения задач только с четными производными Многие прикладные задачи теории пластин, оболочек и конструкций из них математически моделируются с помощью линейных обыкновенных дифференциальных уравнений, содержащих только четные производные. Приведем примеры некоторых из задач, разрешающие уравнения которых содержат только четные производные. Уравнение изгиба балки постоянной жесткости EJ x имеет вид EJ x d 4 y( x) / dx 4  q( x) , где q(x) - поверхностная распределенная нагрузка, у(х) - искомый прогиб балки. Задача о перемещениях мембраны имеет разрешающее уравнение ( 2 / x 2   2 / y 2 )w( x, y)  q / N , где q, N - внешнее поперечное давление и внутреннее нормальное усилие, w(x,y) - искомая функция прогибов. Плоское напряженное состояние пластины, нагруженной контурными силами, лежащими в ее срединной плоскости, описывается уравнением ( 4 / x 4  2 4 / x 2 y 2   4 / y 4 )( x, y)  0 , где  ( x, y ) - функция напряжений Эри. Изгиб изотропных пластин определяется уравнением ( 4 / x 4  2 4 / x 2 y 2   4 / y 4 )w( x, y)  q / D , где q - нормальная к поверхности пластины распределенная внешняя нагрузка, D - цилиндрическая жесткость пластины; w(x,y) искомая функция прогибов. Состояние трехслойной изгибаемой панели, состоящей из тонких несущих слоев и среднего слоя-заполнителя из сот или пенопласта, описывается уравнением  2  2 [1  ( K (1  v) / 2) 2 ]F ( x, y)  q / D,  2 (*)  ( 2 /  2   2 /  2 )(*), D  0.5BH 2 , B  Eh /(1  v 2 ), K  D / C , C  G0 H 48

где D - изгибная жесткость трехслойной панели; E,v,h - модуль Юнга, коэффициент Пуассона и толщина крайних слоев; В - жесткость крайних слоев при растяжении и сжатии; Н - толщина заполнителя, С жесткость заполнителя при поперечном сдвиге; F(x,y) - искомая разрешающая функция. Задача поперечного изгиба ортотропной пластинки сводится к уравнению только с четными производными. В настоящее время для расчетов используются различные варианты уравнений цилиндрических оболочек. Известно, что для круговых цилиндрических оболочек разрешающие уравнения разных авторов Власова, Гольденвейзера, Новожилова, Флюгге, Бейларда, Даревского, Тимошенко - Лява, Доннелла - Власова - Лурье содержат только четные производные. Уравнения в форме Власова содержат только четные производные для изотропных и ортотропных круговых цилиндрических оболочек, нагружаемых силовым воздействием: 1) Полубезмоментная теория, 2) Моментная техническая теория, 3) Общая моментная теория, 4) Общая моментная теория для ортотропной оболочки. Только четные производные содержатся и в разрешающих уравнениях, описывающих деформацию круговых цилиндрических оболочек при воздействии на оболочку температурного поля, постоянного по ее толщине, или поля, характеризующегося перепадом температуры по толщине оболочки: 1) Полубезмоментиая теория, 2) Моментная техническая теория, 3) Общая моментная теория. Дифференциальные уравнения с частными производными сводятся, например, методом Фурье, методом прямых или другими методами разделения переменных к обыкновенным дифференциальным уравнениям, содержащим только четные производные. 10.2. Основы метода Пусть задача описывается следующей системой линейных обыкновенных дифференциальных уравнений в матричном виде y // ( x)  By( x)  f ( x) , 49

где y // ( x)  d 2 y( x) / dx 2 . Здесь B=const, то есть рассмотрим систему обыкновенных дифференциальных уравнений с постоянными коэффициентами. Если в качестве начальных выбираются условия при произвольном значении x  x0 то общее решение задачи Коши (задачи с начальными, а не краевыми условиями) при условии B=const имеет вид [Гантмахер] y( x)  COS ( B1 / 2 ( x  x0 )) y( x0 )  B 1 / 2 SIN( B1 / 2 ( x  x0 )) y / ( x0 ) B 1 / 2 x  SIN(B 1/ 2 ( x  t )) f (t )dt x0 В сокращенных обозначениях: y( x)  C( x  x0 ) y( x0 )  S ( x  x0 ) y / ( x0 )  y0S ( x  x0 ) . Вычисления выполняются по формулам [Гантмахер]: C( x  x0 )  E  (1/ 2!) B( x  x0 ) 2  (1/ 4!) B 2 ( x  x0 ) 4  ... S ( x  x0 )  E( x  x0 )  (1/ 3!) B( x  x0 ) 3  (1/ 5!) B 2 ( x  x0 ) 5  ... x y0S ( x)  B 1 / 2  SIN( B1 / 2 ( x  t )) f (t )dt x0 Из данной формы записи общего решения задачи Коши можно видеть, что для определенных начальных условии конкретное решение задачи зависит от значения двух векторов y( x0 ), y / ( x0 ) . Может быть сформулирована краевая задача. Краевые условия должны быть заданы на значения векторов y(0), y / (0) - на левом крае и y(1), y / (1) - на правом крае. Часто элементы искомой вектор-функции у(х) уравнения не являются искомыми физическими параметрами краевой задачи. В общем случае физические параметры задачи формируются линейными комбинациями элементов вектор-функций y( x), y / ( x) . функции Тогда краевые условия, задаваемые на отдельные физические параметры задачи или на линейные зависимости физических параметров на краях, могут быть представлены в виде: L y(0) l, y / (0) R y(0) r. y / (0) Выражение для вектор-функции дифференцированием выражения для y(x) . 50 y / (x) можно получить

Пользуясь представлением матрицы C( x  x0 ) и матрицы S ( x  x0 ) в виде матричных рядов и выполняя дифференцирование нетрудно показать, что непосредственное их dCOS( B1 / 2 ( x  x0 )) / dx  BB1 / 2 SIN( B1 / 2 ( x  x0 )) , dB 1 / 2 SIN( B1 / 2 ( x  x0 )) / dx  COS ( B1 / 2 ( x  x0 )) . Интеграл x y0S ( x)  B 1 / 2  SIN( B1 / 2 ( x  t )) f (t )dt x0 относится к интегралам, зависящим от параметра, причем с переменными пределами интегрирования и может быть продифференцированным. Собрав полученные производные, запишем здесь искомое выражение для вектор-функции y / ( x)  BS( x  x0 ) y( x0 )  C( x  x0 ) y / ( x0 )  y0C ( x  x0 ) , где y0C ( x  x0 )  B 1 / 2 x  COS (B 1/ 2 ( x  t )) f (t )dt x0 Объединим формулы для y( x0 ), y / ( x0 ) в одну C ( x  x0 ) S ( x  x 0 ) y ( x0 ) y ( x  x0 ) y ( x)   0S / /  BS( x  x0 ) C ( x  x0 ) y ( x0 ) y0C ( x  x0 ) . y ( x) Можем тогда записать C (0  x0 ) S (0  x0 ) y( x0 ) y (0  x0 ) y(0)   0S / /  BS(0  x0 ) C (0  x0 ) y ( x0 ) y0C (0  x0 ) , y (0) C (1  x0 ) S (1  x0 ) y( x0 ) y (1  x0 ) y(1)   0S / /  BS(1  x0 ) C (1  x0 ) y ( x0 ) y0C (1  x0 ) . y (1) Можем записать в совокупности с краевыми условиями по аналогии с методом из главы 3: L C (0  x0 ) S (0  x0 ) y( x0 ) y0S (0  x0 ) y(0)  L (  )l,  BS(0  x0 ) C (0  x0 ) y / ( x0 ) y0C (0  x0 ) y / (0) R C (1  x0 ) S (1  x0 ) y( x0 ) y0S (1  x0 ) y(1)  R (  )r.  BS(1  x0 ) C (1  x0 ) y / ( x0 ) y0C (1  x0 ) y / (1) Тогда получаем одну итоговую систему линейных алгебраических уравнений: 51

L R C (0  x0 ) S (0  x0 )  BS(0  x0 ) C (0  x0 ) C (1  x0 ) S (1  x0 )  BS(1  x0 ) C (1  x0 )  y ( x0 )  y / ( x0 ) lL y0S (0  x0 ) y0C (0  x0 ) rR y0S (1  x0 ) y0C (1  x0 ) Таким образом, записана разрешающая система линейных алгебраических уравнений в смысле и по аналогии с методом из Главы 3, но для дифференциальных уравнений, которые содержат только четные производные. В случае жестких дифференциальных уравнений, содержащих только четные производные, может применяться 1) либо построчное ортонормирование переносимых краевых условий по аналогии с Главой 6, либо 2) может выстраиваться ленточная матрица сопрягаемых участков интервала интегрирования по аналогии с Главой 7; с той разницей, что использовать следует полученные здесь формулы для случая только четных производных в разрешающих дифференциальных уравнениях краевых задач. 52

Приложение 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 все же оказываются необходимыми процедуры ортонормирования. Но на современных 53

персональных компьютерах уже не наблюдаются сплошные скачки значений от больших отрицательных до больших положительных по всему интервалу между краями цилиндра - как это было на компьютерах поколения 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). Компьютеры будут и дальше развиваться такими же темпами как сейчас и это означает, что в самое ближайшее время для подобных расчетов типа расчета неосесимметрично нагруженных оболочек вращения совсем не потребуется применять ортонормирование в рамках предложенного метода «переноса краевых условий», что существенно упрощает программирование метода и увеличивает скорость расчетов не только по сравнению с другими известными методами, но и по сравнению с собственными характеристиками метода «переноса краевых условий» предыдущих лет. 54

ПРОГРАММА НА С++ (РАСЧЕТ ЦИЛИНДРА): //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, где строки матрицы С ортонормированы. 55

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: 56

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); 57

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++){ 58

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]; } 59

} } } } //Вычисление матрицы 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): 60

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 и делается взаимозамена строк 61

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; 62

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};//Соответственно правые части ортонормированного переносимого уравнения с правого края 63

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 64

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 } 65

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; } 66

Приложение 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+(2nn2)/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); 67

} //Задание вектора внешних воздействий в системе ОДУ - вектора 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: 68

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]- 69

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]; 70

} 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; 71

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 первым членом ряда 72

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]; } } } } } //Вычисление 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; } } 73

//Решение СЛАУ размерности 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); 74

} return 0; } int main() { int nn;//Номер гармоники, начиная с 1-й (без нулевой) int nn_last=50;//Номер последней гармоники double Moment[100+1]={0};//Массив физического параметра (момента), что рассчитывается в каждой точке между краями double angle; double start_angle, finish_angle; start_angle=3.14159265359/4; finish_angle=start_angle+(3.14159265359/2); double step=(3.14159265359/2)/100; //step=(3.14159265359/2)/100 - величина шага расчета оболочки - шага интервала интегрирования (должна быть больше нуля, т.е. положительная) double h_div_R;//Величина h/R h_div_R=1.0/200; 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 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 75

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; //Здесь надо первоначально заполнить ненулевыми значениями матрицы и вектора краевых условий U*Y[0]=u_ (слева) и V*Y[100]=v_ (справа) : U[0][0]=nju/tan(start_angle); U[0][1]=1.0; U[0][2]=nju*nn/sin(start_angle); U[0][4]=(1+nju); u_[0]=0.0;//Сила T1 на левом крае равна нулю U[1][0]=-(1-nju)/2/sin(start_angle); U[1][2]=-(1-nju)/2/tan(start_angle); U[1][3]=(1-nju)/2; U[1][4]=-c2*nn*(1-nju)/sin(start_angle)/tan(start_angle); U[1][5]=c2*nn*(1-nju)/sin(start_angle); u_[1]=0.0;//Сила S* на левом краю равна нулю U[2][4]=-nju*nn2/sin(start_angle)/sin(start_angle); U[2][5]=nju/tan(start_angle); U[2][6]=1.0; u_[2]=0;//Момент M1 на левом краю равен нулю U[3][4]=-(3-nju)*nn2/sin(start_angle)/sin(start_angle)/tan(start_angle); U[3][5]=nju+1.0/tan(start_angle)/tan(start_angle)-(nju-2)*nn2/sin(start_angle)/sin(start_angle); U[3][6]=-1.0/tan(start_angle); U[3][7]=-1.0; u_[3]=-sin(nn*gamma)/(nn*gamma);//Сила Q1* на левом крае распределена на угол -gamma +gamma 76

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 angle=start_angle;//начальное значение угловой координаты for(int ii=1;ii<=100;ii++){ angle+=step;//Угловая координата A_perem_coef(nju, c2, nn, angle, A);//Вычисление матрицы А коэффициентов системы ОДУ при данной угловой координате angle exponent(A,(-step),expo_from_minus_step);//Шаг отрицательный (значение шага меньше нуля из-за направления вычисления матричной экспоненты) mat_row_for_partial_vector(A, step, mat_row_for_minus_expo); 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, angle, (-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], 77

//начиная с предпоследней точки - точки с индексом ii=(100-1) используем ii-- (уменьшение индекса точки) angle=finish_angle;//Угловая координата правого края for(int ii=(100-1);ii>=0;ii--){ angle-=step;//Движение справа на лево A_perem_coef(nju, c2, nn, angle, A);//Вычисление матрицы А коэффициентов системы ОДУ при данной угловой координате angle exponent(A,step,expo_from_plus_step);//Шаг положительный (значение шага больше нуля из-за направления вычисления матричной экспоненты) mat_row_for_partial_vector(A, (-step), mat_row_for_plus_expo); 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, angle, 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]; } } //Вычисление момента во всех точках между краями angle=start_angle;//начальное значение угловой координаты for(int ii=0;ii<=100;ii++){ Moment[ii]+=Y[ii][4]*(nju*nn2/sin(angle)/sin(angle))+Y[ii][5]*(nju/tan(angle))+Y[ii][6]*(1.0);//Момент M1 в точке [ii] angle+=step; //U[2][4]=-nju*nn2/sin(start_angle)/sin(start_angle); //U[2][5]=nju/tan(start_angle); //U[2][6]=1.0; Момент } 78

}//ЦИКЛ ПО ГАРМОНИКАМ ЗДЕСЬ ЗАКАНЧИВАЕТСЯ 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; } 79

Приложение 3. Постановка задачи, результаты расчетов и программа на С++ расчета цилиндра – для метода из главы 7 Вычислительные эксперименты проводились в сравнении с методом переноса краевых условий Алексея Юрьевича Виноградова. В этом методе используется построчное ортонормирование. Без ортонормирования в методе переноса краевых условий А.Ю.Виноградова успешно решается задача, например, нагружения цилиндрической оболочки, которая консольно заделана по правому краю и нагружена по левому краю силой, равномерно распределенной по дуге окружности, с отношением длинны к радиусу L/R=2 и с отношением радиуса к толщине R/h=100. Для отношения R/h=200 задача без ортонормирования в методе переноса краевых условий уже не решается, так как выдаются ошибки из-за неустойчивости счета. С применением же ортонормирования в методе переноса краевых условий решаются успешно задачи и для параметров, например, R/h=300, R/h=500, R/h=1000. Новый предлагаемый здесь метод позволяет решать все вышеуказанные тестовые задачи вовсе без применения операций ортонормирования, что значительно упрощает его программирование. Для тестовых расчетов задач с вышеуказанными параметрами новым предлагаемым методом интервал интегрирования разделялся на 10 участков, а между узлами, как и сказано выше, решение находилось как решение задачи Коши. Для решения задач удерживалось 50 гармоник рядов Фурье, так как результат при 50 гармониках уже не отличался от случая удержания 100 гармоник. Скорость же расчета тестовых задач новым предлагаемым методом не меньше, чем методом переноса краевых условий, так как оба метода в тестовых задачах при удержании 50 гармоник рядов Фурье выдавали готовое решение мгновенно после запуска программы на выполнение (на ноутбуке ASUS M51V CPU Duo T5800). В тоже время программирование нового предложенного здесь метода существенно проще, так как нет необходимости программировать процедуры ортонормирования. ПРОГРАММА НА С++ (ЦИЛИНДР): // sopryazhenie.cpp: главный файл проекта. //Решение краевой задачи - цилиндрической оболочки. //Интервал интегрирования разбит на 10 сопрягаемых участков: левый край - точка 0 и правый край точка 10 //БЕЗ ОРТОНОРМИРОВАНИЯ #include "stdafx.h" 80

#include <iostream> #include <conio.h> using namespace std; //Умножение матрицы 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]; } } } //Вычисление матричной экспоненты 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]; 81

} } } } } //Вычисление матрицы 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]; } } } } } 82

//Задание вектора внешних воздействий в системе ОДУ - вектора 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; } } //Решение СЛАУ размерности 88 методом Гаусса с выделением главного элемента int GAUSS(double AA[8*11][8*11], double bb[8*11], double x[8*11]){ double A[8*11][8*11]; double b[8*11]; for(int i=0;i<(8*11);i++){ b[i]=bb[i];//Работать будем с вектором правых частей b, чтобы исходный вектор bb не изменялся при выходе из подпрограммы for(int j=0;j<(8*11);j++){ A[i][j]=AA[i][j];//Работать будем с матрицей А, чтобы исходная матрица АА не менялась при выходе из подпрограммы } } int e;//номер строки, где обнаруживается главный (максимальный) коэффициент в столбце jj double s, t, main;//Вспомогательная величина for(int jj=0;jj<((8*11)-1);jj++){//Цикл по столбцам jj преобразования матрицы А в верхнетреугольную e=-1; s=0.0; main=A[jj][jj]; 83

for(int i=jj;i<(8*11);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*11);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*11);i++){//Приведение к верхнетреугольной матрице for(int j=(jj+1);j<(8*11);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*11)-1]=b[(8*11)-1]/A[(8*11)-1][(8*11)-1];//Первоначальное определение последнего элемента искомого решения х (87-го) for(int i=((8*11)-2);i>=0;i--){//Вычисление елементов решения x[i] от 86-го до 0-го t=0; for(int j=1;j<((8*11)-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.05; //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 84

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 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 Y_many[8*11]={0};// составной вектор из векторов Y(xi) в 11-ти точках с точки 0 (левый край Y(0)) до точки 10 (правый край Y(x10)) double MATRIX_many[8*11][8*11]={0};//матрица СЛАУ double B_many[8*11]={0};// вектор правых частей СЛАУ: MATRIX_many*Y_many=B_many double Y_vspom[8]={0};//вспомогательный вектор double Y_rezult[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; 85

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_ exponent(A,(-step*10),expo_from_minus_step);//Шаг отрицательный (значение шага меньше нуля из-за направления вычисления матричной экспоненты) //x=0.0;//начальное значение координаты - для расчета частного вектора //mat_row_for_partial_vector(A, step, mat_row_for_minus_expo); //Заполнение матрицы коэффициентов СЛАУ MATRIX_many for(int i=0;i<4;i++){ for(int j=0;j<8;j++){ MATRIX_many[i][j]=U[i][j]; MATRIX_many[8*11-4+i][8*11-8+j]=V[i][j]; } B_many[i]=u_[i]; B_many[8*11-4+i]=v_[i]; } for(int kk=0;kk<(11-1);kk++){//(11-1) единичных матриц и матриц EXPO надо записать в MATRIX_many for(int i=0;i<8;i++){ MATRIX_many[i+4+kk*8][i+kk*8]=1.0;//заполнение единичными матрицами for(int j=0;j<8;j++){ MATRIX_many[i+4+kk*8][j+8+kk*8]=expo_from_minus_step[i][j];//заполнение матричными экспонентами } } } //Решение систем линейных алгебраических уравнений GAUSS(MATRIX_many,B_many,Y_many); //Вычисление векторов состояния в 101 точке - левая точка 0 и правая точка 100 exponent(A,step,expo_from_plus_step); 86

for(int i=0;i<11;i++){//заполнение промежуточных точек во всех 10-ти интервалах (всего получим точки от 0 до 100) между 11 узлами for(int j=0;j<8;j++){ Y[0+i*10][j]=Y_many[j+i*8];//в 11-ти узлах векторы беруться из решения СЛАУ - из Y_many } } for(int i=0;i<10;i++){//заполнение промежуточных точек в 10-ти интервалах for(int j=0;j<8;j++){ Y_vspom[j]=Y[0+i*10][j];//начальный вектор для i-го участка, нулевая точка, точка старта i-го участка } mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult); for(int j=0;j<8;j++){ Y[0+i*10+1][j]=Y_rezult[j];//заполнение 1-ой точки интервала Y_vspom[j]=Y_rezult[j];//для следующего шага } mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult); for(int j=0;j<8;j++){ Y[0+i*10+2][j]=Y_rezult[j];//заполнение 2-ой точки интервала Y_vspom[j]=Y_rezult[j];//для следующего шага } mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult); for(int j=0;j<8;j++){ Y[0+i*10+3][j]=Y_rezult[j];//заполнение 3-ой точки интервала Y_vspom[j]=Y_rezult[j];//для следующего шага } mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult); for(int j=0;j<8;j++){ Y[0+i*10+4][j]=Y_rezult[j];//заполнение 4-ой точки интервала Y_vspom[j]=Y_rezult[j];//для следующего шага } mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult); for(int j=0;j<8;j++){ Y[0+i*10+5][j]=Y_rezult[j];//заполнение 5-ой точки интервала Y_vspom[j]=Y_rezult[j];//для следующего шага } mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult); for(int j=0;j<8;j++){ Y[0+i*10+6][j]=Y_rezult[j];//заполнение 6-ой точки интервала Y_vspom[j]=Y_rezult[j];//для следующего шага } mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult); for(int j=0;j<8;j++){ Y[0+i*10+7][j]=Y_rezult[j];//заполнение 7-ой точки интервала Y_vspom[j]=Y_rezult[j];//для следующего шага } mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult); for(int j=0;j<8;j++){ Y[0+i*10+8][j]=Y_rezult[j];//заполнение 8-ой точки интервала 87

Y_vspom[j]=Y_rezult[j];//для следующего шага } mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult); for(int j=0;j<8;j++){ Y[0+i*10+9][j]=Y_rezult[j];//заполнение 9-ой точки интервала Y_vspom[j]=Y_rezult[j];//для следующего шага } } //Вычисление момента во всех точках между краями 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; } 88

Приложение 4. Метод главы 7 и программа для этого метода на С++ на английском языке The simplest method of solution of the stiff boundary value problems – the method of solution of the stiff boundary value problems without orthonormalization offered by Y. I. Vinogradov and A. Y. Vinogradov. Alexei Yurievich Vinogradov (Candidate of science /Physics and Mathematics/) (14 November 2011, Moscow, Russia) The authors: Yuri Ivanovich Vinogradov (1938 year of birth), doctor of science (Physics and Mathematics) and Alexei Yurievich Vinogradov (1970 year of birth), candidate of science (Physics and Mathematics). The idea of overcoming the problems of computing by dividing an integration interval into matching sectors belongs to Y. I. Vinogradov (doctor of science /Physics and Mathematics/). His doctorate thesis was based on that material (including that idea). Proposed realization of this idea (division and matching) through matrix theory formulas i.e. by means of the matrix exponents (or by other words by means of a matrizant – the Cauchy matrix, Cauchy-Krylov functions) belongs to A. Y. Vinogradov (candidate of science /Physics and Mathematics/). AT PRESENT MOMENT THIS METHOD IS THE BEST ONE. The method was finalized in the mid October 2011 and a C++ software program (that is given in this paper as well) was finalized and checked by comparative computations on 14 November 2011. 1. The simplest method of solution of the stiff boundary value problems. Let’s divide a boundary value problem integration interval e.g. for three segments. We will have points (nodes) including the boundaries: x0 , x1 , x2 , x3 . We have the boundary conditions in the following form: UY ( x0 )  u, VY ( x3 )  v. 89

We can write sector matching matrix equations as follows: Y ( x0 )  K ( x0  x1 )Y ( x1 )  Y  ( x0  x1 ) Y ( x1 )  K ( x1  x2 )Y ( x2 )  Y  ( x1  x2 ) Y ( x2 )  K ( x2  x3 )Y ( x3 )  Y  ( x2  x3 ) It is possible to rewrite these formulas in the more convenient form: EY ( x0 )  K ( x0  x1 )Y ( x1 )  Y  ( x0  x1 ) EY ( x1 )  K ( x1  x2 )Y ( x2 )  Y  ( x1  x2 ) EY ( x2 )  K ( x2  x3 )Y ( x3 )  Y  ( x2  x3 ) where E - unit matrix. Then in the combined matrix form we obtain a system of the linear algebraic equations in the following form: U 0 0 0 E  K ( x0  x1 ) 0 0 0 E  K ( x1  x2 ) 0 0 0 E  K ( x 2  x3 ) 0 0 0 V u Y ( x0 )  Y ( x0  x1 ) Y ( x1 )   Y  ( x1  x2 ) Y ( x2 ) Y  ( x 2  x3 ) Y ( x3 ) v This system is solved by the Gauss method by discrimination of the basic element. Solution in the points located between the nodes is obtained by solving the Cauchy problem with initial conditions in the ith node: Y ( x)  K ( x  xi )Y ( xi )  Y  ( x  xi ) It appears to be that there is no need to apply orthonormalization for the boundary value problems for stiff ordinary differential equations. 2. Detailed description of the designations used above. Let’s use as an example a system of differential equations of the cylindrical shell of the rocket – a system of the ordinary differential equations of 8th order (after partial derivatives partitioning by the Fourier method). A system of the linear ordinary differential equations has the following form: Y ( x)  AY ( x)  F ( x) where Y (x) - an unknown vector-function of the 8x1 dimensionality problem, Y (x) - a derivative of the unknown vector-function of 8x1 dimensionality, A - a square matrix of the coefficients of the differential equation of 8x8 dimensionality, F (x ) - a vector-function of the external influence on the system of 8x1 dimensionality. The boundary conditions have the following form: 90

UY (0)  u, VY (1)  v, where Y (0) - is the value of the unknown vector-function at the left boundary x=0 of the 8x1 dimensionality, U - a square horizontal matrix of the coefficients of the boundary conditions of the left boundary of the 4x8 dimensionality, and u - an external influence on the left boundary of the 4x1 dimensionality. Y (1) is the value of the unknown vector-function at the right boundary x=1 of the 8x1 dimensionality, V - a square horizontal matrix of the coefficients of the boundary conditions of the right boundary of the 4x8 dimensionality, and v - an external influence on the left boundary of the 4x1 dimensionality. In case when a system of differential equations has a matrix with constant coefficients A = const the solution of the Cauchy problem has the following form: x Y ( x)  e A( x x 0)Y ( x0 )  e Ax  e  At F (t )dt x0 where e A( x x0)  E  A( x  x0 )  A2 ( x  x0 ) 2 / 2! A3 ( x  x0 ) 3 / 3!..., where E is a unit matrix. Matrix exponent can be named as the Cauchy matrix or matrizant and be denoted in the following form: K ( x  x0 )  K ( x  x0 )  e A( x x0) Then Cauchy problem solution can be written in the following form: Y ( x)  K ( x  x0 )Y ( x0 )  Y  ( x  x0 ) x where Y  ( x  x0 )  e Ax  e  At F (t )dt is a vector of the partial solution of the x0 heterogeneous system of the differential equations. A matrix exponent (the Cauchy matrix) multiplication feature is known from the matrix theory [1]: K ( xi  x0 )  K ( xi  xi 1 )  K ( xi1  xi 2 )  ... K ( x2  x1 )  K ( x1  x0 ) In case when a differential equation system has a matrix with variable coefficients A  A(x) it is offered to find a solution of the Cauchy problem by means of the Cauchy matrix multiplication feature. That is the integration integral is split into small segments and the Cauchy matrixes are computed approximately in these small segments by means of the formula for a constant matrix in the exponent. Then the Cauchy matrixes computed in these small sections are multiplied: K ( xi  x0 )  K ( xi  xi 1 )  K ( xi1  xi 2 )  ... K ( x2  x1 )  K ( x1  x0 ) 91

where the Cauchy matrixes are computed approximately by the following formula: K ( xi 1  xi )  e A( xi) xi  exp( A( xi )  xi ) where xi  xi 1  xi Instead of the formula for computation of the vector of the partial solution of the heterogeneous system of the differential equations in the form [1]: x Y  ( x  x0 )  e Ax  e  At F (t )dt x0 it is proposed to use the following formula for each individual segment of the integration interval: xj Y ( x j  xi )  Y ( x j  xi )  K ( x j  xi )  K ( xi  t ) F (t )dt xi   Correctness of the above mentioned formula is justified by the following: xj Y ( x j  xi )  exp( A( x j  xi ))  exp( A( xi  t ))F (t )dt xi  xj  Y ( x j  xi )   exp( A( x j  xi )) exp( A( xi  t ))F (t )dt xi xj  Y ( x j  xi )   exp( A( x j  xi  xi  t ))F (t )dt xi xj  Y ( x j  xi )   exp( A( x j  t ))F (t )dt xi xj Y ( x j  xi )  exp( Ax j )  exp( At) F (t )dt xi  x Y ( x  xi )  exp( Ax)  exp( At)F (t )dt xi  as was to be justified. Computation of the vector of the partial solution of the heterogeneous system of differential equations is carried out by representing the Cauchy matrix under the integral sign in the form of the series and integrating this series element by element: 92

xj Y ( x j  xi )  Y ( x j  xi )  K ( x j  xi )  K ( xi  t ) F (t )dt  xi xj  K ( x j  xi )  ( E  A( xi  t )  A 2 ( xi  t ) 2 / 2!...)F (t )dt  xi xj xj xj 2  K ( x j  xi )(E  F (t )dt  A  ( xi  t ) F (t )dt  A / 2!  ( xi  t ) 2 F (t )dt  ...). xi xi xi   This formula is correct for a system of differential equations with the constant coefficient matrix A =const. F (t ) vector can be considered in the segment ( x j  xi ) approximately in the form of the constant value F ( xi )  constant that allows its taking off from the integral sign that results in very simple series for computation in the considered segment. In case of the differential equations with variable coefficients the averaged matrix Ai  A( xi ) of the coefficients of the differential equation system can be used in the above mentioned formula for each segment. Let’s consider an option when the integration interval steps are selected as sufficiently small that results in considering F (t ) vector in the segment ( x j  xi ) approximately in the form of the constant value F ( xi )  constantresulting in taking off this vector from the integral signs: xj xj xj Y ( x j  xi )  K ( x j  xi )(E  dt  A  ( xi  t )dt  A / 2!  ( xi  t ) 2 dt  ...)F (t ). xi xi xi  2 It is known that at T=(at+b) we have: T n dt  1 T n 1  const (at n  -1). a(n  1) In our case we have:  (b - t) n dt  1 (b - t) n 1  const (at n  -1). (-1)(n 1) xj Then we obtain  (x xi i  t ) n dt   1 ( xi  x j ) n1 . n 1 Then we obtain the series for computation of the vector of the partial solution of the heterogeneous system of differential equations in the small segment ( x j  xi ) : Y  ( x j  xi )  K ( x j  xi )  ( E  A( xi  x j ) / 2! A2 ( xi  x j ) 2 / 3!...) ( x j  xi )  F ( xi ). 93

If the segment ( x j  xi ) is not a small one it is possible to divide it by subsegments and then it would be possible to propose the following recurrent (iteration) formulas for particular vector computing: We have Y ( x)  K ( x  x0 )Y ( x0 )  Y  ( x  x0 ) . We have the formula for an individual sub-segment as well: xj Y ( x j  xi )  Y ( x j  xi )  K ( x j  xi )  K ( xi  t ) F (t )dt xi   We can write: Y ( x1 )  K ( x1  x0 )Y ( x0 )  Y  ( x1  x0 ) , Y ( x2 )  K ( x2  x1 )Y ( x1 )  Y  ( x2  x1 ) . Let’s substitute Y ( x1 ) in Y ( x2 ) and obtain the following: Y ( x2 )  K ( x2  x1 )[K ( x1  x0 )Y ( x0 )  Y  ( x1  x0 )]  Y  ( x2  x1 )   K ( x2  x1 ) K ( x1  x0 )Y ( x0 )  K ( x2  x1 )Y  ( x1  x0 )  Y  ( x2  x1 ) Let’s compare the obtained expression with the formula: Y ( x2 )  K ( x2  x0 )Y ( x0 )  Y  ( x2  x0 ) and we will obtain evidently that: K ( x2  x0 )  K ( x2  x1 )K ( x1  x0 ) so, we obtain the formula for a particular vector: Y  ( x2  x0 )  K ( x2  x1 )Y  ( x1  x0 )  Y  ( x2  x1 ) Hence, the vectors of the sub-segments Y  ( x1  x0 ),Y  ( x2  x1 ) are not added to each other simply but they are added with involvement of the Cauchy matrix of the sub-segment. Similarly we can write Y ( x3 )  K ( x3  x2 )Y ( x2 )  Y  ( x3  x2 ) and having substituted the formula for Y ( x2 ) in this expression we will obtain: Y ( x3 )  K ( x3  x2 )[K ( x2  x1 ) K ( x1  x0 )Y ( x0 )  K ( x2  x1 )Y  ( x1  x0 )  Y  ( x2  x1 )]   Y  ( x3  x2 )  K ( x3  x2 ) K ( x2  x1 ) K ( x1  x0 )Y ( x0 )   K ( x3  x2 ) K ( x2  x1 )Y  ( x1  x0 )  K ( x3  x2 )Y  ( x2  x1 )  Y  ( x3  x2 ). having compared this obtained expression with the formula: Y ( x3 )  K ( x3  x0 )Y ( x0 )  Y  ( x3  x0 ) we have obtained evidently that: K ( x3  x0 )  K ( x3  x2 )K ( x2  x1 )K ( x1  x0 ) and at the same time we have obtained a formula for a particular vector: Y  ( x3  x0 )  K ( x3  x2 ) K ( x2  x1 )Y  ( x1  x0 )  K ( x3  x2 )Y  ( x2  x1 )  Y  ( x3  x2 ). Hence, a particular vector – a vector of the partial solution of the heterogeneous system of differential equations is computed by this way exactly, for example, the particular vector Y  ( x3  x0 ) is computed similarly in the 94

considered segment ( x3  x0 ) through the computed particular vectors Y  ( x1  x0 ) , Y  ( x2  x1 ) , and Y  ( x3  x2 ) of the corresponding subs-segments ( x1  x0 ) , ( x2  x1 ) and ( x3  x2 ) . Computation accuracy control can be carried out by the following way. For a homogeneous system of differential equations we have: Y ( x)  K ( x  x0 )Y ( x0 ) . We may write the following: Y ( x j )  K ( x j  xi )Y ( xi ) and Y ( xi )  K ( xi  x j )Y ( x j ) Having substituted one formula in another one we obtain: Y ( x j )  K ( x j  xi )Y ( xi )  K ( x j  xi ) K ( xi  x j )Y ( x j ) , i.e. we obtain: Y ( x j )  K ( x j  xi ) K ( xi  x j )Y ( x j ) , however, the latter is possible only in case when K ( x j  xi ) K ( xi  x j )  E - unit matrix, i.e. the matrixes K ( x j  xi ) and K ( xi  x j ) are inverse ones. In other words it has been proved that K 1 ( x j  xi )  K ( xi  x j ) i.e. K 1 ( x j  xi )  K ( xi  x j ) . Hence, computing accuracy can be controlled by means of determination of the accuracy of matrix exponent computation (in other words the Cauchy matrixes or matrizants), that can be checked by having ascertained that the reciprocal inversion of the corresponding matrix exponents is observed in the integration segments: K ( x j  xi ) K ( xi  x j )  E . 3. Computational experiments. Computational experiments have been carried out in comparison with the method of the boundary conditions transfer of Alexey Vinogradov. Line-byline orthonormalization is used in this method. Without use of orthonormalization it is possible by means of the boundary conditions transfer method to solve a problem of cylindrical shell loading that is fixed in cantilever fashion at the right boundary and loaded at the left boundary by the force distributed uniformly by the circular arch with the ratio of the length to the radius L/R=2 and with ratio of the radius to the thickness R/h=100. In case of ratio R/h=200 the problem by means of the 95

method of the boundary conditions transfer without orthonormalozation cannot be solved by this time because there are mistakes resulted in due to counting instability. However, in case of use of orthonormalization in the method of the boundary conditions transfer the problems related to the parameters R/h=300, R/h=500, R/h=1000 can be solved. A new method proposed in this paper allows solving of all above mentioned test problems without use of orthonormalization operation that results in significant simplification of its programming. In case of the test computations of the problems characterized by the above mentioned parameters by means of this new proposed method the integration interval is divided into 10 segments while between the nodes as aforesaid the solution was found as a solution of the Cauchy problem. 50 harmonics of the Fourier series were used for solving the problem since the result in case of usage of 50 harmonics didn’t differ from the case when 100 harmonics were used. Test problem computing speed by means of the proposed method is not less compared to the boundary conditions transfer method because both methods when used for test problems while using 50 harmonics of the Fourier series produced a final solution instantaneously after launching a program (notebook ASUS M51V CPU Duo T5800). At the same time programming of this newly proposed method is significantly simpler because there is no need in orthonormalization procedure programming. LIST OF REFERENCES 1. Gantmaher F.R. Matrix theory. – M.: Nauka, 1988. – 548 p. C++ program // sopryazhenie.cpp: main file of the project. //Solution of the boundary value problem – a cylindrical shell problem. //Integration interval is divided into 10 matching segments: left boundary – point 0 and right boundary – point 10. //WITHOUT ORTHONORMALIZATION #include "stdafx.h" #include <iostream> #include <conio.h> using namespace std; //Multiplication of A matrix by b vector and obtain 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]; } } 96

} //Computation of the matrix exponent EXP=exp(A*delta_x) void exponent(double A[8][8], double delta_x, double EXP[8][8]) { //n – number of the terms of the series in the exponent, m – a counter of the number of the terms of the series (m<=n) int n=100, m; double E[8][8]={0}, TMP1[8][8], TMP2[8][8]; int i,j,k; //E – unit matrix – the first term of the series of the exponent 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; //initial filling-in of the auxiliary array TMP1 – the previous term of the series for follow-up multiplication //and initial filling-in of the exponent by the first term of the series for(i=0;i<8;i++) { for(j=0;j<8;j++) { TMP1[i][j]=E[i][j]; EXP[i][j]=E[i][j]; } } //series of EXP exponent computation starting from the 2nd term of the series (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;//taken out beyond the cycle of multiplication of the row by the column TMP2[i][j]/=(m-1);// taken out beyond the cycle of multiplication of the row by the column EXP[i][j]+=TMP2[i][j]; } } //filling-in of the auxiliary array TMP1 for computing the next term of the series TMP2 in the next step of the cycle by m if (m<n) { for(i=0;i<8;i++) { for(j=0;j<8;j++) { TMP1[i][j]=TMP2[i][j]; } } } } } //computation of the matrix MAT_ROW in the form of the matrix series for follow-up use //when computing a vector - partial vector – a vector of the partial solution of the heterogeneous system of the ordinary differential equations at the step delta x void mat_row_for_partial_vector(double A[8][8], double delta_x, double MAT_ROW[8][8]) { 97

//n – number of the terms of the series in MAT_ROW, m – a counter of the number of the terms of the series (m<=n) int n=100, m; double E[8][8]={0}, TMP1[8][8], TMP2[8][8]; int i,j,k; //E – unit matrix – the first term of the series 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; //initial filling-in of the auxiliary array TMP1 – the previous term of the series for following multiplication //and initial filling-in of MAT_ROW by the first term of the series 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]; } } //a series of computation of MAT_ROW starting from the second term of the series (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]; } } //filling-in of the auxiliary array TMP1 for computing the next term of the series – TMP2 in the next step of the cycle by m if (m<n) { for(i=0;i<8;i++) { for(j=0;j<8;j++) { TMP1[i][j]=TMP2[i][j]; } } } } } //specifying the external influence vector in the system of ordinary differential equations – POWER vector: 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; } //computation of the vector – ZERO (particular case) vector of the partial solution //heterogeneous system of differential equations in the segment of interest: 98

void partial_vector(double vector[8]){ for(int i=0;i<8;i++){ vector[i]=0.0; } } //computation of the vector – the vector of the partial solution of the heterogeneous system of differential equations in the segment of interest 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};//External influence vector on the shell double REZ[8]={0}; double REZ_2[8]={0}; power_vector_for_partial_vector(x_, POWER_);//Computing POWER_ at coordinate x_ mat_on_vect(mat_row, POWER_, REZ);//Multiplication of the matrix mat_row by POWER vector and obtain REZ vector mat_on_vect(expo_, REZ, REZ_2);// Multiplication of matrix expo_ by vector REZ and obtain vector REZ_2 for(int i=0;i<8;i++){ vector[i]=REZ_2[i]*delta_x; } } //Solution of SLAE of 88 dimensionality by the Gauss method with discrimination of the basic element int GAUSS(double AA[8*11][8*11], double bb[8*11], double x[8*11]){ double A[8*11][8*11]; double b[8*11]; for(int i=0;i<(8*11);i++){ b[i]=bb[i];//we will operate with the vector of the и right parts to provide that initial vector bb would not change when exiting the subprogram for(int j=0;j<(8*11);j++){ A[i][j]=AA[i][j];//we will operate with A matrix to provide that initial AA matrix would not change when exiting the subprogram } } int e;//number of the row where main (maximal) coefficient in the column jj is found double s, t, main;//Ancillary quantity for(int jj=0;jj<((8*11)-1);jj++){//Cycle by columns jj of transformation of A matrix into upper-triangle one e=-1; s=0.0; main=A[jj][jj]; for(int i=jj;i<(8*11);i++){//there is a number of у row where main (maximal) element is placed in the column jj and row interchanging is made if ((A[i][jj]*A[i][jj])>s) {//Instead of multiplication (potential minus sign is deleted) it could be possible to use a function by abs() module e=i; s=A[i][jj]*A[i][jj]; } } if (e<0) { cout<<"Mistake "<<jj<<"\n"; return 0; } if (e>jj) {//If the main element isn’t placed in the row with jj number but is placed in the row with у number main=A[e][jj]; for(int j=0;j<(8*11);j++){//interchanging of two rows with e and jj numbers t=A[jj][j]; A[jj][j]=A[e][j]; A[e][j]=t; 99

} t=b[jj]; b[jj]=b[e]; b[e]=t; } for(int i=(jj+1);i<(8*11);i++){//reduction to the upper-triangle matrix for(int j=(jj+1);j<(8*11);j++){ A[i][j]=A[i][j]-(1/main)*A[jj][j]*A[i][jj];//re-calculation of the coefficients of the row i>(jj+1) } b[i]=b[i]-(1/main)*b[jj]*A[i][jj]; A[i][jj]=0.0;//nullified elements of the row under diagonal element of A matrix } }//Cycle by jj columns of transformation of A matrix into upper-triangle one x[(8*11)-1]=b[(8*11)-1]/A[(8*11)-1][(8*11)-1];//initial determination of the last element of the desired solution x (87th) for(int i=((8*11)-2);i>=0;i--){//Computation of the elements of the solution x[i] from 86th to 0th t=0; for(int j=1;j<((8*11)-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;//Number of the harmonic starting from the 1st (without zero one) int nn_last=50;//Number of the last harmonic double Moment[100+1]={0};//An array of the physical parameter (momentum) that is computed in each point between the boundaries double step=0.05; //step=(L/R)/100 – step size of shell computation – a step of integration interval (it should be over zero, i.e. be positive) double h_div_R;//Value of h/R h_div_R=1.0/100; double c2; c2=h_div_R*h_div_R/12;//Value of h*h/R/R/12 double nju; nju=0.3; double gamma; gamma=3.14159265359/4;//The force distribution angle by the left boundary //printing to files: 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++){ //A CYCLE BY HARMONICS STARTING FROM THE 1 st HARMONIC (EXCEPT ZERO ONE) 100

double x=0.0;//A coordinate from the left boundary – it is needed in case of heterogeneous system of the ODE for computing the particular vector FF double expo_from_minus_step[8][8]={0};//The matrix for placement of the exponent in it at the step of (0-x1) type double expo_from_plus_step[8][8]={0};// The matrix for placement of the exponent in it in the step of (x1-0) type double mat_row_for_minus_expo[8][8]={0};//the auxiliary matrix for particular vector computing when moving at step of (0-x1) type double mat_row_for_plus_expo[8][8]={0};// the auxiliary matrix for particular vector computing when moving at step of (x1-0) type double U[4][8]={0};//The matrix of the boundary conditions of the left boundary of the dimensionality 4x8 double u_[4]={0};//Dimensionality 4 vector of the external influence for the boundary conditions of the left boundary double V[4][8]={0};//The boundary conditions matrix of the right boundary of the dimensionality 4x8 double v_[4]={0};// The dimensionality 4 vector of the external influence for the boundary conditions for the right boundary double Y[100+1][8]={0};//The array of the vector-solutions of the corresponding linear algebraic equations system (in each point of the interval between the boundaries): MATRIXS*Y=VECTORS double A[8][8]={0};//Matrix of the coefficients of the system of ODE double FF[8]={0};//Vector of the particular solution of the heterogeneous ODE at the integration interval sector double Y_many[8*11]={0};// a composite vector consisting of the vectors Y(xi) in 11 points from point 0 (left boundary Y(0) to the point 10 (right boundary Y(x10)) double MATRIX_many[8*11][8*11]={0};//The matrix of the system of the linear algebraic equations double B_many[8*11]={0};// a vector of the right parts of the SLAE: MATRIX_many*Y_many=B_many double Y_vspom[8]={0};//an auxiliary vector double Y_rezult[8]={0};//an auxiliary vector double nn2,nn3,nn4,nn5,nn6,nn7,nn8;//Number of the nn harmonic raised to corresponding powers nn2=nn*nn; nn3=nn2*nn; nn4=nn2*nn2; nn5=nn4*nn; nn6=nn4*nn2; nn7=nn6*nn; nn8=nn4*nn4; //Filling-in of non-zero elements of the A matrix of the coefficients of ODE system 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; //Here firstly it is necessary to make filling-in by non-zero values of the matrix and boundary conditions vector U*Y[0]=u_ (on the left) and V*Y[100]=v_ (on the right): U[0][1]=1.0; U[0][2]=nn*nju; U[0][4]=nju; u_[0]=0.0;//Force T1 at the left boundary is equal to zero U[1][0]=-(1-nju)/2*nn; U[1][3]=(1-nju)/2; U[1][5]=(1-nju)*nn*c2; u_[1]=0.0;//Force S* at the left boundary is equal to zero U[2][4]=-nju*nn2; U[2][6]=1.0; u_[2]=0;//Momentum M1 at the left boundary is equal to zero U[3][5]=(2-nju)*nn2; U[3][7]=-1.0; u_[3]=-sin(nn*gamma)/(nn*gamma);//Force Q1* at the left boundary is distributed at the angle -gamma +gamma V[0][0]=1.0; v_[0]=0.0;// The right boundary displacement u is equal to zero V[1][2]=1.0; v_[1]=0.0;//The right boundary displacement v is equal to zero V[2][4]=1.0; v_[2]=0.0;//The right boundary displacement w is equal to zero V[3][5]=1.0; v_[3]=0.0;//The right boundary rotation angle is equal to zero //Here initial filling-in of U*Y[0]=u_ и V*Y[100]=v_ terminates 101

exponent(A,(-step*10),expo_from_minus_step);//A negative step (step value is less than zero due to direction of matrix exponent computation) //x=0.0;//the initial value of coordinate – for partial vector computation //mat_row_for_partial_vector(A, step, mat_row_for_minus_expo); //Filling-in of the SLAE coefficients matrix MATRIX_many for(int i=0;i<4;i++){ for(int j=0;j<8;j++){ MATRIX_many[i][j]=U[i][j]; MATRIX_many[8*11-4+i][8*11-8+j]=V[i][j]; } B_many[i]=u_[i]; B_many[8*11-4+i]=v_[i]; } for(int kk=0;kk<(11-1);kk++){//(11-1) unit matrix and EXPO matrix should be written into MATRIX_many for(int i=0;i<8;i++){ MATRIX_many[i+4+kk*8][i+kk*8]=1.0;//filling-in by unit matrixes for(int j=0;j<8;j++){ MATRIX_many[i+4+kk*8][j+8+kk*8]=-expo_from_minus_step[i][j];//fillingin by matrix exponents } } } //Solution of the system of linear algebraic equations GAUSS(MATRIX_many,B_many,Y_many); //Computation of the state vectors in 101 points – left point 0 and right point 100 exponent(A,step,expo_from_plus_step); for(int i=0;i<11;i++){//passing points filling-in in all 10 intervals (we will obtain points from 0 to 100) between 11 nodes for(int j=0;j<8;j++){ Y[0+i*10][j]=Y_many[j+i*8];//in 11 nodes the vectors are taken from SLAE solutions – from Y many } } for(int i=0;i<10;i++){//passing points filling-in in 10 intervals for(int j=0;j<8;j++){ Y_vspom[j]=Y[0+i*10][j];//the initial vector for ith segment, zero point, starting point of the ith segment } mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult); for(int j=0;j<8;j++){ Y[0+i*10+1][j]=Y_rezult[j];//the first point of the interval filling-in Y_vspom[j]=Y_rezult[j];//for the next step } mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult); for(int j=0;j<8;j++){ Y[0+i*10+2][j]=Y_rezult[j];//filling-in of the 2nd point of the interval Y_vspom[j]=Y_rezult[j];//for the next step } mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult); 102

for(int j=0;j<8;j++){ Y[0+i*10+3][j]=Y_rezult[j];//filling-in of the 3rd point of the interval Y_vspom[j]=Y_rezult[j];//for the next step } mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult); for(int j=0;j<8;j++){ Y[0+i*10+4][j]=Y_rezult[j];//filing-in of the 4th point of the interval Y_vspom[j]=Y_rezult[j];//for the next step } mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult); for(int j=0;j<8;j++){ Y[0+i*10+5][j]=Y_rezult[j];// filing-in of the 5th point of the interval Y_vspom[j]=Y_rezult[j];// for the next step } mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult); for(int j=0;j<8;j++){ Y[0+i*10+6][j]=Y_rezult[j];// filing-in of the 6th point of the interval Y_vspom[j]=Y_rezult[j];// for the next step } mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult); for(int j=0;j<8;j++){ Y[0+i*10+7][j]=Y_rezult[j];// filing-in of the 7th point of the interval Y_vspom[j]=Y_rezult[j];// for the next step } mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult); for(int j=0;j<8;j++){ Y[0+i*10+8][j]=Y_rezult[j];// filing-in of the 8th point of the interval Y_vspom[j]=Y_rezult[j];// for the next step } mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult); for(int j=0;j<8;j++){ Y[0+i*10+9][j]=Y_rezult[j];// filing-in of the 9th point of the interval Y_vspom[j]=Y_rezult[j];// for the next step } } //Computation of the momentum in all points between the boundaries for(int ii=0;ii<=100;ii++){ Moment[ii]+=Y[ii][4]*(-nju*nn2)+Y[ii][6]*1.0;//Momentum M1 in the point [ii] //U[2][4]=-nju*nn2; U[2][6]=1.0; u_[2]=0;//Momentum M1 } }//CYCLE BY HARMONICS TERMINATES HERE 103

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; } My web-site related to stiff boundary value problems solution methods: www.vinogradov-design.narod.ru/math.html 104

Приложение 5. Графическое предложение метода численного решения дифференциальных уравнений Эти графические эскизы были сформулированы А.Ю. Виноградовым до окончания МГТУ им. Н.Э.Баумана – до 1993 года. После чего были выложены в интернет в апреле 2007 года. Предположим, что дано дифференциальное уравнение y/(x)=a(x)y(x)+f(x) с начальными условиями y(0)=y0. Или дано аналогичное матричное дифференциальное уравнение с начальными условиями. Пусть решение надо найти не аналитически, а численно. Есть множество методов разных авторов как это можно сделать численно. Простейший и самый ошибочный метод таков: из дифференциального уравнения на основе начальных условий высчитывается значение производной в начальной точке y/(0)=a(0)y(0)+f(0), а потом делается шаг до x1 в предположении, что искомая функция на участке от 0 до x1 моделируется отрезком прямой линии и так далее шаг за шагом. Геометрически это выглядит вот так, как известно: Разные авторы разных численных методов использовали разные приемы получения формул численного интегрирования с идеей уменьшения ошибочности результатов вычислений. Далее предлогаются простейшие эскизы с геометрическими идеями уменьшения ошибочности результатов численного интегрирования, что, конечно же, может быть сформулировано не только графически, но и формулами для выполнения приближенных вычислений: 105

106

Список опубликованных работ 1. Виноградов А.Ю. Вычисление начальных векторов для численного решения краевых задач // - Деп. в ВИНИТИ, 1994. -N2073В94. -15 с. 2. Виноградов А.Ю., Виноградов Ю.И. Совершенствование метода прогонки Годунова для задач строительной механики // Изв. РАН Механика твердого тела, 1994. -№4. -С. 187-191. 3. Виноградов А.Ю. Вычисление начальных векторов для численного решения краевых задач // Журнал вычислительной математики и математической физики, 1995. -Т.З5. -№1. -С. 156-159. 4. Виноградов А.Ю. Численное моделирование произвольных краевых условий для задач строительной механики тонкостенных конструкций // Тез. докладов Белорусского Конгресса по теоретической и прикладной механике "Механика-95". Минск, 6-11 февраля 1995 г., Гомель: Изд- во ИММС АНБ, 1995. -С63-64. 5. Vinogradov A.Yu. Numerical modeling of boundary conditions in deformation problems of structured material in thin wall constructions // International Symposium "Advances in Structured and Heterogeneous Continua II". Book Absfracts. August, 14-16, 1995, Moscow, Russia. -P.51. 6. Виноградов А.Ю. Численное моделирование краевых условий в задачах деформирования тонкостенных конструкций из композиционных материалов // Механика Композиционных Материалов и Конструкций, 1995. -T.I. -N2. - С. 139-143. 7. Виноградов А.Ю. Приведение краевых задач механики элементов приборных устройств к задачам Коши для выбранной точки // Прикладная механика в приборных устройствах. Меж вуз. сб. научных трудов. -Москва: МИРЭА, 1996. 8. Виноградов А.Ю. Модификация метода Годунова// Труды Международной научно-технической конференции "Современные проблемы машиноведения", Гомель: ГПИ им. П.О. Сухого, 1996. - С.39-41. 9. Виноградов А.Ю., Виноградов Ю.И. Методы переноса сложных краевых условий для жестких дифференциальных уравнений строительной механики // Труды Международной конференции «Ракетно-космическая техника: фундаментальные проблемы механики и теплообмена», Москва, 1998. 10. Виноградов А.Ю. Метод решения краевых задач путем переноса условий с краев интервала интегрирования в произвольную 107

точку // Тез. докладов Международной конференции "Актуальные проблемы механики оболочек", Казань, 2000.-С. 176. 11. Виноградов Ю.И., Виноградов А.Ю., Гусев Ю.А. Метод переноса краевых условий для дифференциальных уравнений теории оболочек // Труды Международной конференции "Актуальные проблемы механики оболочек", Казань, 2000. -С. 128-132. 12. Виноградов А.Ю., Виноградов Ю.И. Метод переноса краевых условий функциями Коши-Крылова для жестких линейных обыкновенных дифференциальных уравнений. // ДАН РФ, – М.: 2000, т. 373, №4, с. 474-476. 13. Виноградов А.Ю., Виноградов Ю.И. Функции Коши-Крылова и алгоритмы решения краевых задач теории оболочек // ДАН РФ, - М.: 2000. -Т.375.-№3.-С. 331-333. 14. Виноградов А.Ю. Численные методы переноса краевых условий // Журнал "Математическое моделирование", изд-во РАН, Институт математического моделирования, - М.: 2000, Т. 12 , № 7, с. 3-6. 15. Виноградов А.Ю., Гусев Ю.А. Перенос краевых условий функциями Коши-Крылова в задачах строительной механики // Тез. докладов VIII Всероссийского съезда по теоретической и прикладной механике, Пермь, 2001. -С. 109-110. 16. Виноградов Ю.И., Виноградов А.Ю. Перенос краевых условий функциями Коши-Крылова // Тез. докладов Международной конференции "Dynamical System Modeling and Stability Investigation""DSMSI-2001", Киев, 2001. 17. Виноградов А.Ю., Виноградов Ю.И., Гусев Ю.А, Клюев Ю.И. Перенос краевых условий функциями Коши-Крылова и его свойства // Изв. РАН МТТ, №2. 2001. -С.155-161. 18. Виноградов Ю.И., Виноградов А.Ю., Гусев Ю.А.Численный метод переноса краевых условий для жестких дифференциальных уравнений строительной механики // Журнал "Математическое моделирование", изд-во РАН, Институт математического моделирования, - М.: 2002, Т. 14, №9, с.3-8. 19. Виноградов Ю.И., Виноградов А.Ю. Простейший метод решения жестких краевых задач // Фундаментальные исследования. – 2014. – № 12–12. – С. 2569-2574; 20. Виноградов Ю.И., Виноградов А.Ю. Решение жестких краевых задач строительной механики (расчет оболочек составных и со шпангоутами) методом Виноградовых (без ортонормирования) // Современные проблемы науки и образования. – 2015. - №1; 108

21. Виноградов Ю.И., Виноградов А.Ю. Уточненный метод Виноградовых переноса краевых условий в произвольную точку интервала интегрирования для решения жестких краевых задач // Современные наукоемкие технологии. – 2016. – № 11-2. – С. 226-232; 22. Виноградов А.Ю. Методы решения жестких и нежестких краевых задач: монография // Волгоград: Изд-во ВолГУ, 2016. – 128 с. 23. Виноградов А.Ю. Программа (код) на С++ решения жесткой краевой задачи методом А.Ю. Виноградова // exponenta.ru/educat/systemat/vinogradov/index.asp 24. Виноградов А.Ю. Методы А.Ю.Виноградова решения краевых задач, в том числе жестких краевых задач // exponenta.ru/educat/systemat/vinogradov/index2.asp 25. Виноградов А.Ю. Расчет составных оболочек и оболочек со шпангоутами простейшим методом решения жестких краевых задач (без ортонормирования) // exponenta.ru/educat/systemat/vinogradov/index3.asp 26. Виноградов А.Ю. Методы решения жестких и нежестких краевых задач – народная докторская // exponenta.ru/educat/systemat/vinogradov/index5.asp 27. Виноградов А.Ю. Методы решения жестких и нежестких краевых задач (9 методов) // exponenta.ru/educat/systemat/vinogradov/index6.asp 109

Научное издание Алексей Юрьевич Виноградов Численные методы решения жестких и нежестких краевых задач Публикуется в авторской редакции Рассылка обязательных экземпляров осуществлена в соответствии с действующими нормативными актами Подписано в печать 05.01.2017г. Формат 60х84/16. Усл. п.л. 2,6. Тираж 1000. Заказ 204. National Research Москва, Россия 2017

Рецензии:

Авторизуйтесь, чтобы добавить рецензию

- у работы пока нет рецензий -

Отзывы:

Авторизуйтесь, чтобы оставить отзыв