Вычисление вектора частного решения неоднородной системы дифференциальных уравнений. Методика "переноса краевых условий" в произвольную точку интервала интегрирования. Расчет обратной матрицы. Замена метода численного интегрирования Рунге-Кутта.
Аннотация к работе
На примере системы дифференциальных уравнений цилиндрической оболочки ракеты - системы обыкновенных дифференциальных уравнений 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 )
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 с.