Программа на С++ решения жесткой краевой задачи методом А.Ю.Виноградова.

1. Введение.
            Метод проверен и его эффективность выше, чем эффективность (скорость счета) метода ортогональной прогонки С.К.Годунова - до 2 порядков (в 100 раз). Для разных типов задач преимущество по скорости разное. А для некоторых типов задач метод А.Ю.Виноградова не требует применять ортонормирование вовсе.
Рассмотрим метод на примере системы дифференциальных уравнений цилиндрической оболочки ракеты – системы обыкновенных дифференциальных уравнений 8-го порядка (после разделения частных производных).
Система линейных обыкновенных дифференциальных уравнений имеет вид:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru],
где ** – искомая вектор-функция задачи размерности 8х1, ** – производная искомой вектор-функции размерности 8х1, ** – квадратная матрица коэффициентов дифференциального уравнения размерности 8х8, *** – вектор-функция внешнего воздействия на систему размерности 8х1.
Здесь и далее вектора обозначаем жирным шрифтом вместо черточек над буквами
Краевые условия имеют вид:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru]
где      ** – значение искомой вектор-функции на левом крае х=0 размерности 8х1, ** – прямоугольная горизонтальная матрица коэффициентов краевых условий левого края размерности 4х8, * – вектор внешних воздействий на левый край размерности 4х1,
** – значение искомой вектор-функции на правом крае х=1 размерности 8х1, **– прямоугольная горизонтальная матрица коэффициентов краевых условий правого края размерности 4х8, * – вектор внешних воздействий на правый край размерности 4х1.
В случае, когда система дифференциальных уравнений имеет матрицу с постоянными коэффициентами **=const, решение задачи Коши имеет вид [1]:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru],
 
где [здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru], где - это единичная матрица.
Матричная экспонента ещё может называться матрицей Коши или матрициантом и может обозначаться в виде:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru].
Тогда решение задачи Коши может быть записано в виде:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru],
где [здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru] это вектор частного решения неоднородной системы дифференциальных уравнений.
2. Случай переменных коэффициентов.
Из теории матриц [1] известно свойство перемножаемости матричных экспонент (матриц Коши):
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru]
 
В случае, когда система дифференциальных уравнений имеет матрицу с переменными коэффициентами , решение задачи Коши предлагается искать при помощи свойства перемножаемости матриц Коши. То есть интервал интегрирования разбивается на малые участки и на малых участках матрицы Коши приближенно вычисляются по формуле для постоянной матрицы в экспоненте. А затем матрицы Коши, вычисленные на малых участках, перемножаются:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru],
где матрицы Коши приближенно вычисляются по формуле:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru], где ***.
3. Вычисление вектора частного решения неоднородной системы дифференциальных уравнений.
Вместо формулы для вычисления вектора частного решения неоднородной системы дифференциальных уравнений в виде [1]:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru]
предлагается использовать следующую формулу для каждого отдельного участка интервала интегрирования:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru].
Правильность приведенной формулы подтверждается следующим:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru],
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru],
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru],
***,
***,
***,
что и требовалось подтвердить.
Вычисление вектора частного решения неоднородной системы дифференциальных уравнений производиться при помощи представления матрицы Коши под знаком интеграла в виде ряда и интегрирования этого ряда поэлементно:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru]
Эта формула справедлива для случая системы дифференциальных уравнений с постоянной матрицей коэффициентов =const.
Вектор ** может рассматриваться на участке ** приближенно в виде постоянной величины **, что позволяет вынести его из под знака интеграла, что приводит к совсем простому ряду для вычислений на рассматриваемом участке.
Для случая дифференциальных уравнений с переменными коэффициентами в приведенной выше формуле для каждого участка может использоваться осредненная матрица ** коэффициентов системы дифференциальных уравнений.
Рассмотрим вариант, когда шаги интервала интегрирования выбираются достаточно малыми, что позволяет рассматривать вектор ** на участке ** приближенно в виде постоянной величины **, что позволяет вынести этот вектор из под знаков интегралов:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru]
Известно, что при T=(at+b) имеем [здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru]
В нашем случае имеем [здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru]
Тогда получаем [здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru].
Тогда получаем ряд для вычисления вектора частного решения неоднородной системы дифференциальных уравнений на малом участке **:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru]
Приведем формулы вычисления вектора частного решения, например, *** на рассматриваемом участке *** через вектора частного решения ***, ***, *** соответствующих подучастков ***, ***, ***.
Имеем [здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru].
Также имеем формулу для отдельного подучастка:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru].
Можем записать:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru],
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru].
Подставим *** в *** и получим:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru]
***.
Сравним полученное выражение с формулой:
***
и получим, очевидно, что:
***
и для частного вектора получаем формулу:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru].
То есть вектора подучастков [здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru] не просто складываются друг с другом, а с участием матрицы Коши подучастка.
            Аналогично запишем *** и подставим сюда формулу для *** и получим:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru]
            Сравнив полученное выражение с формулой:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru]
очевидно, получаем, что:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru]
и вместе с этим получаем формулу для частного вектора:
***
То есть именно так и вычисляется частный вектор – вектор частного решения неоднородной системы дифференциальных уравнений, то есть так вычисляется, например, частный вектор *** на рассматриваемом участке *** через вычисленные частные вектора ***, ***, *** соответствующих подучастков ***, ***, ***.
4. Метод «переноса краевых условий» в произвольную точку интервала интегрирования.
Полное решение системы дифференциальных уравнений имеет вид
***.
Или можно записать:
***.
Подставляем это выражение для *** в краевые условия левого края и получаем:
***,
***,
***.
Или получаем краевые условия, перенесенные в точку ***:
***,
где     ***     и      ***.
Далее запишем аналогично
***
И подставим это выражение для ***в перенесенные краевые условия точки ***:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru],
***,
***.
Или получаем краевые условия, перенесенные в точку ***:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru],
где     ***     и      ***.
И так в точку *** переносим матричное краевое условие с левого края и таким же образом переносим матричное краевое условие с правого края.
Покажем шаги переноса краевых условий с правого края.
Можем записать:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru]
Подставляем это выражение для *** в краевые условия правого края и получаем:
***,
***,
***
Или получаем краевые условия правого края, перенесенные в точку ***:
***,
где     ***     и      ***.
Далее запишем аналогично
***
И подставим это выражение для ***в перенесенные краевые условия точки ***:
***,
***,
***.
Или получаем краевые условия, перенесенные в точку ***:
***,
где     ***     и      ***.
И так во внутреннюю точку *** интервала интегрирования переносим матричное краевое условие, как показано, и с левого края и таким же образом переносим матричное краевое условие с правого края и получаем:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru],
***.
Из этих двух матричных уравнений с прямоугольными горизонтальными матрицами коэффициентов очевидно получаем одну систему линейных алгебраических уравнений с квадратной матрицей коэффициентов:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru].
5. Случай «жестких» дифференциальных уравнений.
В случае «жестких» дифференциальных уравнений предлагается применять построчное ортонормирование матричных краевых условий в процессе их переноса в рассматриваемую точку. Для этого формулы ортонормирования систем линейных алгебраических уравнений можно взять в [2].
То есть, получив
***
применяем к этой группе линейных алгебраических уравнений построчное ортонормирование и получаем эквивалентное матричное краевое условие:
***.
 
И теперь уже в это проортонормированное построчно уравнение подставляем
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru].
И получаем
***,
***.
Или получаем краевые условия, перенесенные в точку ***:
***,
где     ***[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru]     и      ***.
Теперь уже к этой группе линейных алгебраических уравнений применяем построчное ортонормирование и получаем эквивалентное матричное краевое условие:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru]
И так далее.
И аналогично поступаем с промежуточными матричными краевыми условиями, переносимыми с правого края в рассматриваемую точку.
В итоге получаем систему линейных алгебраических уравнений с квадратной матрицей коэффициентов, состоящую из двух независимо друг от друга поэтапно проортонормированных матричных краевых условий. Эта система решается методом Гаусса с выделением главного элемента для получения решения **** в рассматриваемой точке ***:
***.
6. Контроль точности вычислений.
Для однородной системы дифференциальных уравнений имеем:
***.
Можем записать:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru] и
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru]
Подставляем одну формулу в другую и получаем:
**,
то есть получаем:
***,
но последнее возможно только когда
*** - единичная матрица,
то есть матрицы *** и *** взаимообратны.
То есть доказано, что
***,
то есть
***.
            Таким образом, получаем, что точность вычислений можно контролировать при помощи определения точности вычисления матричных экспонент (или иначе говоря - матриц Коши или матрициантов), что можно проверять, удостоверяясь, что на участках интегрирования выполняется условие взаимной обратности соответствующих матричных экспонент:
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru].
7. Применяемые формулы ортонормирования.
Взято из [2]. Пусть дана система линейных алгебраических уравнений порядка n:
**=*.
 
Здесь над векторами поставим черточки вместо их обозначения жирным шрифтом.
Будем рассматривать строки матрицы * системы как векторы:
*=(*,*,…,*).
 
Ортонормируем эту систему векторов.
Первое уравнение системы **=* делим на *.
При этом получим:
**+**+…+**=*,    *=(*,*,…,*),
где[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru]=, =1.
 
Второе уравнение системы заменяется на:[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru]
8. Вычислительные эксперименты.
В качестве проверочных задач использовалась схема консольно закрепленных цилиндрической и сферической оболочек с параметрами R/h=50, 100, 200. Длина цилиндрической оболочки рассматривалась L/R=2, а угловые координаты сферической оболочки рассматривались от p/4 до 3p/4. На свободном крае рассматривалось нормальное к поверхности оболочек погонное усилие, равномерно распределенное в интервале [-p/4, p/4]. В качестве среды программирования использовалась система MicrosoftVisualStudio 2010 (VisualC++).
Первоначально метод был предложен и обсчитывался в кандидатской диссертации А.Ю.Виноградова в 1993-1995 годах. Тогда оказалось, что без использования ортонормирования в рамках метода успешно решаются задачи осесимметрично нагруженных оболочек вращения. Расчеты тогда выполнялись на компьютере поколения 286. Задачи же неосесимметричного нагружения оболочек вращения можно было решать на компьютерах поколения 286 только с применением процедур построчного ортонормирования - как это и предлагалось в рамках метода. Без процедур ортонормирования в неосесимметричных случаях выдавались только ошибочные графики, представлявшие собой хаотично скачущие большие отрицательные и большие положительные значения, например, изгибающего обезразмеренного момента М1.
Современные компьютеры имеют значительно более совершенное внутреннее устройство и более точные внутренние операции с числами, чем это было в 1993-1995 годах. Поэтому было интересно рассмотреть возможность расчета неосесимметрично нагруженных оболочек, например, цилиндров, на современном аппаратном и программном обеспечении в рамках предложенного метода «переноса краевых условий» совсем без использования процедур построчного ортонормирования.
Оказалось, что неосесимметрично нагруженные цилиндры при некоторых параметрах на современных компьютерах уже можно решать в рамках предложенного метода «переноса краевых условий» совсем без применения операций построчного ортонормирования. Это, например, при параметрах цилиндра L/R=2 и R/h=100.
При параметрах цилиндра L/R=2 и R/h=200 все же оказываются необходимыми процедуры ортонормирования. Но на современных персональных компьютерах уже не наблюдаются сплошные скачки значений от больших отрицательных до больших положительных по всему интервалу между краями цилиндра - как это было на компьютерах поколения 286. В частном случае L/R=2 и R/h=200 наблюдаются лишь незначительные скачки в районе максимума изгибающего обезразмеренного момента М1 на левом крае и небольшой скачек обезразмеренного момента М1 на правом крае.
Приводятся графики изгибающего обезразмеренного момента М1:
- слева приводятся графики, полученные при использовании операций построчного ортонормирования на каждом из 100 шагов, на которые разделялся участок между краями цилиндра,
- справа приводятся графики, полученные совсем без применения операций построчного ортонормирования.
Следует сказать, что в качестве расчетной среды использовалась 32-х битная операционная система WindowsXP и среда программирования MicrosoftVisualStudio 2010 (VisualC++) использовалась в тех же рамках 32-х битной организации операций с числами. Параметры компьютера такие: ноутбук ASUSM51V (CPUDuoT5800).
Компьютеры будут и дальше развиваться такими же темпами как сейчас и это означает, что в самое ближайшее время для подобных расчетов типа расчета неосесимметрично нагруженных оболочек вращения совсем не потребуется применять ортонормирование в рамках предложенного метода «переноса краевых условий», что существенно упрощает программирование метода и увеличивает скорость расчетов не только по сравнению с другими известными методами, но и по сравнению с собственными характеристиками метода «переноса краевых условий» предыдущих лет.
[здесь идёт математическая формула - чтобы её увидеть нужно скать архив с данной работой, это абслютно бесплатно на referatwork.ru]

СПИСОК ЛИТЕРАТУРЫ

1. Гантмахер Ф.Р. Теория матриц. – М.: Наука, 1988. – 548 с.
2. Березин И.С., Жидков Н.П. Методы вычислений, том II, Государственное издательство физико-математической литературы, Москва, 1962 г., 635 с.

P.S. Метод Алексея Юрьевича Виноградова «переноса краевых условий» первоначально был опубликован в Межвузовском сборнике МИРЭА (кажется в 1995 году). МИРЭА это Московский институт радиотехники, электроники и автоматики. Точное название и год выхода стать можно посмотреть в Ленинской библиотеке в списке литературы моей диссертации. Там у меня только одна статья в МИРЭА. К сожалению, на руках у меня нет экземпляра моей кандидатской диссертации, поэтому не могу привести точное название статьи, но называется она, кажется, что-то вроде «Метод приведения краевых задач к задаче Коши».