Исправленные методы А.Ю. Виноградова: решения краевых задач, в том числе жестких краевых задач - Научная работа

бесплатно 0
4.5 172
Вычисление вектора частного решения неоднородной системы дифференциальных уравнений. Методика "переноса краевых условий" в произвольную точку интервала интегрирования. Расчет обратной матрицы. Замена метода численного интегрирования Рунге-Кутта.


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

Система линейных обыкновенных дифференциальных уравнений имеет вид: Y (x) = A(x) • Y(x) F(x), где Y(x) - искомая вектор-функция задачи размерности 8х1, Y (x) - производная искомой вектор-функции размерности 8х1, A(x) - квадратная матрица коэффициентов дифференциального уравнения размерности 8х8, F(x) - вектор-функция внешнего воздействия на систему размерности 8х1.

Здесь и далее вектора обозначаем жирным шрифтом вместо черточек над буквами

Краевые условия имеют вид: U•Y(0) = u, V•Y(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, решение задачи Коши имеет вид [Гантмахер]: Y(x) = e • Y(x ) e • e • F(t) dt, где e = E A(x-x ) A (x-x ) /2! A (x-x ) /3! …, где E это единичная матрица.

Матричная экспонента еще может называться матрицей Коши или матрициантом и может обозначаться в виде: K(x<x ) = K(x - x ) = e .

Тогда решение задачи Коши может быть записано в виде: Y(x) = K(x<x ) • Y(x ) Y*(x<x ) , где

Y*(x<x ) = e • e • F(t) dt это вектор частного решения неоднородной системы дифференциальных уравнений

1. Случай переменных коэффициентов

Этот вариант рассмотрения переменных коэффициентов проверялся в кандидатской диссертации.

Из теории матриц [Гантмахер] известно свойство перемножаемости матричных экспонент (матриц Коши): e = e • e • … • e • e , K(x <x ) = K(x <x ) • K(x <x ) • … • K(x <x ) • K(x <x ).

В случае, когда система дифференциальных уравнений имеет матрицу с переменными коэффициентами A=A(x), решение задачи Коши предлагается искать при помощи свойства перемножаемости матриц Коши. То есть интервал интегрирования разбивается на малые участки и на малых участках матрицы Коши приближенно вычисляются по формуле для постоянной матрицы в экспоненте. А затем матрицы Коши, вычисленные на малых участках, перемножаются: K(x <x ) = K(x <x ) • K(x <x ) • … • K(x <x ) • K(x <x ), где матрицы Коши приближенно вычисляются по формуле: K(x <x ) = e , где ?x = x - x .

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

Эта очень простая формула еще не обсчитана на компьютерах. Вместо нее обсчитывалась (а может быть и не обсчитывалась, не знаю) значительно ранее выведенная (выведенная моим отцом) и гораздо более сложная формула, приведенная в: Численный метод переноса краевых условий для жестких дифференциальных уравнений строительной механики Журнал "ММ", Том: 14 (2002), Номер: 9, 3 стр. 1409-003r.pdf

Вместо формулы для вычисления вектора частного решения неоднородной системы дифференциальных уравнений в виде [Гантмахер]: Y*(x<x ) = e • e • F(t) dt предлагается использовать следующую формулу для каждого отдельного участка интервала интегрирования: Y*(x <x ) = Y*(x - x ) = K(x - x ) • K(x - t) • F(t) dt .

Правильность приведенной формулы подтверждается следующим: Y*(x - x ) = e • e • F(t) dt , Y*(x - x ) = e •e • F(t) dt , Y*(x - x ) = e • F(t) dt , Y*(x - x ) = e • F(t) dt , Y*(x - x ) = e • e • F(t) dt , Y*(x<x ) = e • e • F(t) dt, что и требовалось подтвердить.

Вычисление вектора частного решения системы дифференциальных уравнений производиться при помощи представления матрицы Коши под знаком интеграла в виде ряда и интегрирования этого ряда поэлементно: Y*(x <x ) = Y*(x - x ) = K(x - x ) • K(x - t) • F(t) dt =

= K(x - x ) • (E A(x - t) A (x - t) /2! … ) • F(t) dt =

= K(x - x ) • (E F(t) dt A• (x - t) • F(t) dt A /2! • (x - t) • F(t) dt … ) .

Эта формула справедлива для случая системы дифференциальных уравнений с постоянной матрицей коэффициентов A=const.

Вектор F(t) может рассматриваться на участке (x - x ) приближенно в виде постоянной величины F(х )=constant, что позволят вынести его из под знака интеграла, что приводит к совсем простому ряду для вычислений на рассматриваемом участке.

Для случая дифференциальных уравнений с переменными коэффициентами в приведенной выше формуле для каждого участка может использоваться осредненная матрица А: A =А(х ) коэффициентов системы дифференциальных уравнений.

Приведем (итерационные или рекуррентные) формулы вычисления вектора частного решения, например, Y*(x <x ) на рассматриваемом участке (x <x ) через вектора частного решения Y*(x <x ), Y

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

Система линейных обыкновенных дифференциальных уравнений имеет вид: Y (x) = A(x) • Y(x) F(x), где Y(x) - искомая вектор-функция задачи размерности 8х1, Y (x) - производная искомой вектор-функции размерности 8х1, A(x) - квадратная матрица коэффициентов дифференциального уравнения размерности 8х8, F(x) - вектор-функция внешнего воздействия на систему размерности 8х1.

Здесь и далее вектора обозначаем жирным шрифтом вместо черточек над буквами

Краевые условия имеют вид: U•Y(0) = u, V•Y(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, решение задачи Коши имеет вид [Гантмахер]: Y(x) = e • Y(x ) e • e • F(t) dt, где e = E A(x-x ) A (x-x ) /2! A (x-x ) /3! …, где E это единичная матрица.

Матричная экспонента еще может называться матрицей Коши или матрициантом и может обозначаться в виде: K(x<x ) = K(x - x ) = e .

Тогда решение задачи Коши может быть записано в виде: Y(x) = K(x<x ) • Y(x ) Y*(x<x ) , где

Y*(x<x ) = e • e • F(t) dt это вектор частного решения неоднородной системы дифференциальных уравнений

1. Случай переменных коэффициентов

Этот вариант рассмотрения переменных коэффициентов проверялся в кандидатской диссертации.

Из теории матриц [Гантмахер] известно свойство перемножаемости матричных экспонент (матриц Коши): e = e • e • … • e • e , K(x <x ) = K(x <x ) • K(x <x ) • … • K(x <x ) • K(x <x ).

В случае, когда система дифференциальных уравнений имеет матрицу с переменными коэффициентами A=A(x), решение задачи Коши предлагается искать при помощи свойства перемножаемости матриц Коши. То есть интервал интегрирования разбивается на малые участки и на малых участках матрицы Коши приближенно вычисляются по формуле для постоянной матрицы в экспоненте. А затем матрицы Коши, вычисленные на малых участках, перемножаются: K(x <x ) = K(x <x ) • K(x <x ) • … • K(x <x ) • K(x <x ), где матрицы Коши приближенно вычисляются по формуле: K(x <x ) = e , где ?x = x - x .

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

Эта очень простая формула еще не обсчитана на компьютерах. Вместо нее обсчитывалась (а может быть и не обсчитывалась, не знаю) значительно ранее выведенная (выведенная моим отцом) и гораздо более сложная формула, приведенная в: Численный метод переноса краевых условий для жестких дифференциальных уравнений строительной механики Журнал "ММ", Том: 14 (2002), Номер: 9, 3 стр. 1409-003r.pdf

Вместо формулы для вычисления вектора частного решения неоднородной системы дифференциальных уравнений в виде [Гантмахер]: Y*(x<x ) = e • e • F(t) dt предлагается использовать следующую формулу для каждого отдельного участка интервала интегрирования: Y*(x <x ) = Y*(x - x ) = K(x - x ) • K(x - t) • F(t) dt .

Правильность приведенной формулы подтверждается следующим: Y*(x - x ) = e • e • F(t) dt , Y*(x - x ) = e •e • F(t) dt , Y*(x - x ) = e • F(t) dt , Y*(x - x ) = e • F(t) dt , Y*(x - x ) = e • e • F(t) dt , Y*(x<x ) = e • e • F(t) dt, что и требовалось подтвердить.

Вычисление вектора частного решения системы дифференциальных уравнений производиться при помощи представления матрицы Коши под знаком интеграла в виде ряда и интегрирования этого ряда поэлементно: Y*(x <x ) = Y*(x - x ) = K(x - x ) • K(x - t) • F(t) dt =

= K(x - x ) • (E A(x - t) A (x - t) /2! … ) • F(t) dt =

= K(x - x ) • (E F(t) dt A• (x - t) • F(t) dt A /2! • (x - t) • F(t) dt … ) .

Эта формула справедлива для случая системы дифференциальных уравнений с постоянной матрицей коэффициентов A=const.

Вектор F(t) может рассматриваться на участке (x - x ) приближенно в виде постоянной величины F(х )=constant, что позволят вынести его из под знака интеграла, что приводит к совсем простому ряду для вычислений на рассматриваемом участке.

Для случая дифференциальных уравнений с переменными коэффициентами в приведенной выше формуле для каждого участка может использоваться осредненная матрица А: A =А(х ) коэффициентов системы дифференциальных уравнений.

Приведем (итерационные или рекуррентные) формулы вычисления вектора частного решения, например, Y*(x <x ) на рассматриваемом участке (x <x ) через вектора частного решения Y*(x <x ), Y*(x <x ), Y*(x <x ) соответствующих подучастков (x <x ), (x <x ), (x <x ).

Имеем: Y(x) = K(x<x ) • Y(x ) Y*(x<x ) , Также имеем формулу для отдельного подучасточка: Y*(x <x ) = Y*(x - x ) = K(x - x ) • K(x - t) • F(t) dt.

Можем записать: Y(x ) = K(x <x ) • Y(x ) Y*(x <x ) , Y(x ) = K(x <x ) • Y(x ) Y*(x <x ) .

Подставим Y(x ) в Y(x ) и получим: Y(x ) = K(x <x ) • [K(x <x ) • Y(x ) Y*(x <x )] Y*(x <x ) =

= K(x <x ) • K(x <x ) • Y(x ) K(x <x ) • Y*(x <x ) Y*(x <x ).

Сравним полученное выражение с формулой: Y(x ) = K(x <x ) • Y(x ) Y*(x <x ) и получим, очевидно, что K(x <x ) = K(x <x ) • K(x <x ) и, самое главное здесь - для частного вектора получаем формулу: Y*(x <x ) = K(x <x ) • Y*(x <x ) Y*(x <x ).

То есть вектора подучастков Y*(x <x ) и Y*(x <x ) не просто складываются друг с другом, а с участием матриц Коши подучастков.

Аналогично запишем: Y(x ) = K(x <x ) • Y(x ) Y*(x <x )

И подставим сюда формулу для Y(x ) и получим: Y(x ) = K(x <x ) • [K(x <x ) • K(x <x ) • Y(x ) K(x <x ) • Y*(x <x ) Y*(x <x )] Y*(x <x ) = K(x <x ) • K(x <x ) • K(x <x ) • Y(x ) K(x <x ) • K(x <x ) • Y*(x <x ) K(x <x ) • Y*(x <x ) Y*(x <x ).

Сравнив полученное выражение с формулой: Y(x ) = K(x <x ) • Y(x ) Y*(x <x )

K(x <x ) = K(x <x ) • K(x <x ) • K(x <x ) и вместе с этим получаем формулу для частного вектора: Y*(x <x ) = K(x <x ) • K(x <x ) • Y*(x <x ) K(x <x ) • Y*(x <x ) Y*(x <x ).

То есть именно так (по своеобразным рекуррентным формулам) и вычисляется частный вектор - вектор частного решения неоднородной системы дифференциальных уравнений, то есть так вычисляется, например, частный вектор Y*(x <x ) всего участка (x <x ) на основе вычисленных векторов Y*(x <x ), Y*(x <x ), Y*(x <x ) подучастков (x <x ), (x <x ), (x <x ).

3. Метод "переноса краевых условий" в произвольную точку интервала интегрирования

Метод обсчитан на компьютерах. По нему уже сделано 3 кандидатских физмат диссертации.

Метод подходит для любых краевых задач. А для "жестких" краевых задач показано, что метод считает быстрее, чем метод С.К.Годунова до 2-х порядков (в 100 раз), а для некоторых "жестких" краевых задач не требует ортонормирования вовсе. Численный метод переноса краевых условий для жестких дифференциальных уравнений строительной механики Журнал "ММ", Том: 14 (2002), Номер: 9, 3 стр. 1409-003r.pdf

Полное решение системы дифференциальных уравнений имеет вид

Y(x) = K(x<x ) • Y(x ) Y*(x<x ) .

Или можно записать: Y(0) = K(0<x ) • Y(x ) Y*(0<x ) .

Подставляем это выражение для Y(0) в краевые условия левого края и получаем: U•Y(0) = u, U•[ K(0<x ) • Y(x ) Y*(0<x ) ] = u, [ U• K(0<x ) ] • Y(x ) = u - U•Y*(0<x ) .

Или получаем краевые условия, перенесенные в точку x : U • Y(x ) = u , где

U = [ U• K(0<x ) ] и u = u - U•Y*(0<x ) .

Далее запишем аналогично

Y(x ) = K(x <x ) • Y(x ) Y*(x <x )

И подставим это выражение для Y(x ) в перенесенные краевые условия точки x

U • Y(x ) = u , U • [ K(x <x ) • Y(x ) Y*(x <x ) ] = u , [ U • K(x <x ) ] • Y(x ) = u - U • Y*(x <x ) , Или получаем краевые условия, перенесенные в точку x : U • Y(x ) = u , где

U = [ U • K(x <x ) ] и u = u - U • Y*(x <x ) .

Покажем перенос краевых условий с правого края.

Можно записать: Y(1) = K(1<x ) • Y(x ) Y*(1<x ) .

Подставляем это выражение для Y(1) в краевые условия правого края и получаем: V•Y(1) = v, V•[ K(1<x ) • Y(x ) Y*(1<x ) ] = v, [ V• K(1<x )] • Y(x ) = v - V• Y*(1<x ).

Или получаем краевые условия, перенесенные в точку x : V • Y(x ) = v , где

V • = [V• K(1<x )] и v = v - V• Y*(1<x ).

Далее запишем аналогично

Y(x ) = K(x <x ) • Y(x ) Y*(x <x )

И подставим это выражение для Y(x ) в перенесенные краевые условия точки x

V • Y(x ) = v , V • [K(x <x ) • Y(x ) Y*(x <x ) ] = v , [V • K(x <x )] • Y(x ) = v - V • Y*(x <x ).

Или получаем краевые условия, перенесенные в точку x : V • Y(x ) = v , где

V • = [V • K(x <x )] и v = v - V • Y*(x <x ).

И так в точку x переносим матричное краевое условие с левого края и таким же образом переносим матричное краевое условие с правого края и получаем: U • Y(x ) = u , V • Y(x ) = v .

Из этих двух матричных уравнений с прямоугольными горизонтальными матрицами коэффициентов очевидно получаем одну систему линейных алгебраических уравнений с квадратной матрицей коэффициентов: • Y(x ) = .

А в случае "жестких" дифференциальных уравнений предлагается применять построчное ортонормирование матричных краевых условий в процессе их переноса в рассматриваемую точку. Для этого формулы ортонормирования систем линейных алгебраических уравнений можно взять в [Березин, Жидков]. То есть, получив

U • Y(x ) = u , применяем к этой группе линейных алгебраических уравнений построчное ортонормирование и получаем эквивалентное матричное краевое условие:

U • Y(x ) = u .

И теперь уже в это проортонормированное построчно уравнение подставляем

Y(x ) = K(x <x ) • Y(x ) Y*(x <x ) .

И получаем

U • [ K(x <x ) • Y(x ) Y*(x <x ) ] = u , [ U • K(x <x ) ] • Y(x ) = u - U • Y*(x <x ) , Или получаем краевые условия, перенесенные в точку x : U • Y(x ) = u , где

U = [ U • K(x <x ) ] и u = u - U • Y*(x <x ) .

Теперь уже к этой группе линейных алгебраических уравнений применяем построчное ортонормирование и получаем эквивалентное матричное краевое условие: U • Y(x ) = u .

И так далее.

И аналогично поступаем с промежуточными матричными краевыми условиями, переносимыми с правого края в рассматриваемую точку.

В итоге получаем систему линейных алгебраических уравнений с квадратной матрицей коэффициентов, состоящую из двух независимо друг от друга поэтапно проортонормированных матричных краевых условий. Эта система решается методом Гаусса с выделением главного элемента для получения решения Y(x ) в рассматриваемой точке x : • Y(x ) = .

4. Программа на С расчета цилиндрической оболочки

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

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

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

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

Список литературы
1. Гантмахер Ф.Р. Теория матриц. - М.: Наука, 1988. - 548 с.

2. Березин И.С., Жидков Н.П. Методы вычислений, том II, Государственное издательство физико-математической литературы, Москва, 1962 г., 635 с.

3. www.vinogradov-design.narod.ru/math.html

4. www.vinogradov-best.narod.ru

5. www.ALEXEIVINOGRADOV.narod.ru

6. www.VINOGRADOVALEXEI.narod.ru

7. www.Vinogradov-Alexei.narod.ru

8. www.Vinogradov-math.narod.ru

Размещено на .ru
Заказать написание новой работы



Дисциплины научных работ



Хотите, перезвоним вам?