--PAGE_BREAK--4 Метод «переноса краевых условий» в произвольную точку интервала интегрирования
Метод обсчитан на компьютерах. По нему уже сделано 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) .
И так в точку 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) = .
5 Второй вариант метода «переноса краевых условий» в произвольную точку интервала интегрирования
Этот вариант метода еще не обсчитан на компьютерах.
Предложено выполнять интегрирование по формулам теории матриц [Гантмахер] сразу от некоторой внутренней точки интервала интегрирования к краям:
Y(0) = K(0←x) ∙ Y(x) + Y*(0←x) ,
Y(1) = K(1←x) ∙ Y(x) + Y*(1←x) .
Подставим эти формулы в краевые условия и получим:
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) .
и
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:
[ U∙ K(0←x) ] ∙ Y(x) = u— U∙Y*(0←x) ,
[ V∙ K(1←x) ] ∙ Y(x) = v— V∙Y*(1←x) .
Эти уравнения аналогично объединяются в одну систему линейных алгебраических уравнений с квадратной матрицей коэффициентов для нахождения решения Y(x) в любой рассматриваемой точке x:
∙ Y(x) = .
В случае «жестких» дифференциальных уравнений предлагается следующий алгоритм.
Используем свойство перемножаемости матриц Коши:
K(x←x) = K(x←x) ∙ K(x←x) ∙ … ∙ K(x←x) ∙ K(x←x)
и запишем выражения для матриц Коши, например, в виде:
K(0←x) = K(0←x) ∙ K(x←x) ∙ K(x←x),
K(1←x) = K(1←x) ∙ K(x←x) ∙ K(x←x) ∙ K(x←x),
Тогда перенесенные краевые условия можно записать в виде:
[ U∙ K(0←x) ∙ K(x←x) ∙ K(x←x) ] ∙ Y(x) = u— U∙Y*(0←x) ,
[ V∙ K(1←x) ∙ K(x←x) ∙ K(x←x) ∙ K(x←x) ] ∙ Y(x) = v— V∙Y*(1←x)
или в виде:
[ U∙ K(0←x) ∙ K(x←x) ∙ K(x←x) ] ∙ Y(x) = u* ,
[ V∙ K(1←x) ∙ K(x←x) ∙ K(x←x) ∙ K(x←x) ] ∙ Y(x) = v* .
Тогда рассмотрим левое перенесенное краевое условие:
[ U∙ K(0←x) ∙ K(x←x) ∙ K(x←x) ] ∙ Y(x) = u* ,
[ U∙ K(0←x) ] ∙ { K(x←x) ∙ K(x←x) ∙ Y(x) } = u* ,
[ матрица ] ∙ { вектор } = вектор .
Эту группу линейных алгебраических уравнений можно подвергнуть построчному ортонормированию, которое сделает строчки [матрицы] ортонормированными, {вектор} затронут не будет, а вектор получит преобразование. То есть получим:
[ U∙ K(0←x) ] ∙ { K(x←x) ∙ K(x←x) ∙ Y(x) } = u* .
Далее последовательно можно записать:
[[ U∙ K(0←x) ] ∙ K(x←x) ] ∙ { K(x←x) ∙ Y(x) } = u* ,
[ матрица ] ∙ { вектор } = вектор .
Аналогично и эту группу линейных алгебраических уравнений можно подвергнуть построчному ортонормированию, которое сделает строчки [матрицы] ортонормированными, {вектор} затронут не будет, а вектор получит преобразование. То есть получим:
[[ U∙ K(0←x) ] ∙ K(x←x) ] ∙ { K(x←x) ∙ Y(x) } = u* ,
Далее аналогично можно записать:
[[[ U∙ K(0←x) ] ∙ K(x←x) ] ∙ K(x←x) ] ∙ { Y(x) } = u* ,
[ матрица ] ∙ { вектор} = вектор .
Аналогично и эту группу линейных алгебраических уравнений можно подвергнуть построчному ортонормированию, которое сделает строчки [матрицы] ортонормированными, {вектор} затронут не будет, а вектор получит преобразование. То есть получим:
[[[ U∙ K(0←x) ]∙ K(x←x) ] ∙ K(x←x) ] ∙ Y(x) = u* .
Аналогично можно проортонормировать матричное уравнение краевых условий и для правого края независимо от левого края.
Далее проортонормированные уравнения краевых условий:
[ U∙ K(0←x) ] ∙ Y(x) = u* ,
[ V∙ K(1←x) ] ∙ Y(x) = v*
как и ранее объединяются в одну обычную систему линейных алгебраических уравнений с квадратной матрицей коэффициентов для нахождения искомого вектора Y(x) :
∙ Y(x) = .
продолжение
--PAGE_BREAK--6 Метод дополнительных краевых условий
Этот метод еще не обсчитан на компьютерах.
Запишем на левом крае ещё одно уравнение краевых условий:
M∙ Y(0) = m.
В качестве строк матрицы M можно взять те краевые условия, то есть выражения тех физических параметров, которые не входят в параметры краевых условий левого края L или линейно независимы с ними. Это вполне возможно, так как у краевых задач столько независимых физических параметров какова размерность задачи, а в параметры краевых условий входит только половина физических параметров задачи. То есть, например, если рассматривается задача об оболочке ракеты, то на левом крае могут быть заданы 4 перемещения. Тогда для матрицы М можно взять параметры сил и моментов, которых тоже 4, так как полная размерность такой задачи – 8. Вектор m правой части неизвестен и его надо найти и тогда можно считать, что краевая задача решена, то есть сведена к задаче Коши, то есть найден вектор Y(0) из выражения:
∙ Y(0) = ,
то есть вектор Y(0) находится из решения системы линейных алгебраических уравнений с квадратной невырожденной матрицей коэффициентов, состоящей из блоков U и M.
Аналогично запишем на правом крае ещё одно уравнение краевых условий:
N∙ Y(0) = n,
где матрица N записывается из тех же соображений дополнительных линейно независимых параметров на правом крае, а вектор nнеизвестен.
Для правого края тоже справедлива соответствующая система уравнений:
∙ Y(1) = .
Запишем Y(1) = K(1←0) ∙Y(0) + Y*(1←0) и подставим в последнюю систему линейных алгебраических уравнений:
∙ [ K(1←0) ∙Y(0) + Y*(1←0) ] = ,
∙ K(1←0) ∙Y(0) = — ∙ Y*(1←0),
∙ K(1←0) ∙Y(0) = ,
∙ K(1←0) ∙Y(0) = .
Запишем вектор Y(0) через обратную матрицу:
Y(0) = ∙
и подставим в предыдущую формулу:
∙ K(1←0) ∙ ∙ = .
Таким образом, мы получили систему уравнений вида:
В ∙ = ,
где матрица В известна, векторы u и sизвестны, а векторы mи tнеизвестны.
Разобьем матрицу В на естественные для нашего случая 4 блока и получим:
∙ = ,
откуда можем записать, что
В11 ∙ u+ B12 ∙ m= s,
B21 ∙ u+ B22 ∙ m= t.
Следовательно, искомый вектор m вычисляется по формуле:
m= B12∙ (s – B11∙ u).
А искомый вектор n вычисляется через вектор t:
t = B21 ∙ u+ B22 ∙ m,
n= t+ N∙ Y*(1←0).
В случае «жестких» дифференциальных уравнений предлагается выполнять поочередное построчное ортонормирование.
Запишем приведенную выше формулу
∙ K(1←0) ∙ ∙ =
в виде:
∙ K(1←x2) ∙ K(x2←x1) ∙ K(x1←0) ∙ ∙ = .
Эту формулу можно записать в виде разделения левой части на произведение матрицы на вектор:
[ ∙ K(1←x2) ] ∙ { K(x2←x1) ∙ K(x1←0) ∙ ∙ } =
[ матрица ] ∙ { вектор } = вектор
Эту группу линейных алгебраических уравнений можно подвергнуть построчному ортонормированию, которое сделает строчки [матрицы] ортонормированными, {вектор} затронут не будет, а вектор получит преобразование. То есть получим:
[ ∙ K(1←x2) ] { K(x2←x1) ∙ K(x1←0) ∙ ∙ } =
Здесь следует сказать, что подвектор t подвергать преобразованию не нужно, так как невозможно, так как его первоначальное значение не известно. Но подвектор t нам оказывается и не нужен для решения задачи.
Далее запишем:
[[ ∙ K(1←x2) ] ∙ K(x2←x1)] { K(x1←0) ∙ ∙ } =
[ матрица ] { вектор } = вектор
Аналогично и эту группу линейных алгебраических уравнений можно подвергнуть построчному ортонормированию, которое сделает строчки [матрицы] ортонормированными, {вектор} затронут не будет, а вектор получит преобразование. То есть получим:
[[ ∙ K(1←x2) ] K(x2←x1)] { K(x1←0) ∙ } = .
И так далее.
В результате поочередного ортонормирования получим:
В ∙ = ,
∙ = .
Следовательно, искомый вектор m вычисляется по формуле:
m= B12 ∙ (s – B11∙ u).
7 Формула для начала счета методом прогонки С.К.Годунова
Эта формула обсчитана на компьютерах в кандидатской диссертации.
Рассмотрим проблему метода прогонки С.К.Годунова.
редположим, что рассматривается оболочка ракеты. Это тонкостенная труба. Тогда система линейных обыкновенных дифференциальных уравнений будет 8-го порядка, матрица A(x) коэффициентов будет иметь размерность 8х8, искомая вектор-функция Y(x) будет иметь размерность 8х1, а матрицы краевых условий будут прямоугольными горизонтальными размерности 4х8.
Тогда в методе прогонки С.К.Годунова для такой задачи решение ищется в следующем виде:
Y(x) = Y(x) c + Y(x) c + Y(x) c + Y(x) c + Y*(x),
или можно записать в матричном виде:
Y(x) = Y(x) ∙ c + Y*(x),
где векторы Y(x), Y(x), Y(x), Y(x) – это линейно независимые вектора-решения однородной системы дифференциальных уравнений, а вектор Y*(x) – это вектор частного решения неоднородной системы дифференциальных уравнений.
Здесь Y(x)=|| Y(x), Y(x), Y(x), Y(x) || это матрица размерности 8х4, а c это соответствующий вектор размерности 4х1из искомых констант c,c,c,c.
Но вообще то решение для такой краевой задачи с размерностью 8 (вне рамок метода прогонки С.К.Годунова) может состоять не из 4 линейно независимых векторов Y(x), а полностью из всех 8 линейно независимых векторов-решений однородной системы дифференциальных уравнений:
Y(x)=Y(x)c+Y(x)c+Y(x)c+Y(x)c+
+Y(x)c+Y(x)c+Y(x)c+Y(x)c+Y*(x),
И как раз трудность и проблема метода прогонки С.К.Годунова и состоит в том, что решение ищется только с половиной возможных векторов и констант и проблема в том, что такое решение с половиной констант должно удовлетворять условиям на левом крае (стартовом для прогонки) при всех возможных значениях констант, чтобы потом найти эти константы из условий на правом крае.
То есть в методе прогонки С.К.Годунова есть проблема нахождения таких начальных значений Y(0), Y(0), Y(0), Y(0), Y*(0) векторов Y(x), Y(x), Y(x), Y(x), Y*(x), чтобы можно было начать прогонку с левого края x=0, то есть чтобы удовлетворялись условия U∙Y(0) = uна левом крае при любых значениях констант c,c,c,c.
Обычно эта трудность «преодолевается» тем, что дифференциальные уравнения записываются не через функционалы, а через физические параметры и рассматриваются самые простейшие условия на простейшие физические параметры, чтобы начальные значения Y(0), Y(0), Y(0), Y(0), Y*(0) можно было угадать. То есть задачи со сложными краевыми условиями так решать нельзя: например, задачи с упругими условиями на краях.
Ниже предлагается формула для начала вычислений методом прогонки С.К.Годунова.
Выполним построчное ортонормирование матричного уравнения краевых условий на левом крае:
U∙Y(0) = u,
где матрица U прямоугольная и горизонтальная размерности 4х8.
В результате получим эквивалентное уравнение краевых условий на левом крае, но уже с прямоугольной горизонтальной матрицей U размерности 4х8, у которой будут 4 ортонормированные строки:
U∙Y(0) = u,
где в результате ортонормирования вектор u преобразован в вектор u.
Как выполнять построчное ортонормирование систем линейных алгебраических уравнений можно посмотреть в [Березин, Жидков].
Дополним прямоугольную горизонтальную матрицу U до квадратной невырожденной матрицы W:
W= ,
где матрица М размерности 4х8 должна достраивать матрицу U до невырожденной квадратной матрицы W размерности 8х8.
В качестве строк матрицы М можно взять те краевые условия, то есть выражения тех физических параметров, которые не входят в параметры левого края или линейно независимы с ними. Это вполне возможно, так как у краевых задач столько независимых физических параметров какова размерность задачи, то есть в данном случае их 8 штук и если 4 заданы на левом крае, то ещё 4 можно взять с правого края.
Завершим ортонормирование построенной матрицы W, то есть выполним построчное ортонормирование и получим матрицу W размерности 8х8 с ортонормированными строками:
W = .
Можем записать, что
Y(0) = (М)транспонированная = М.
Тогда, подставив в формулу метода прогонки С.К.Годунова, получим:
Y(0) = Y(0) ∙с + Y*(0)
или
Y(0) = М∙с + Y*(0).
Подставим эту последнюю формулу в краевые условия U∙Y(0) = u и получим:
U∙ [ М∙с + Y*(0) ]= u.
Отсюда получаем, что на левом крае константы c уже не на что не влияют, так как
U∙ М = 0 и остается только найти Y*(0) из выражения:
U∙ Y*(0) = u.
Но матрица U имеет размерность 4х8 и её надо дополнить до квадратной невырожденной, чтобы найти вектор Y*(0) из решения соответствующей системы линейных алгебраических уравнений:
∙ Y*(0) = ,
где 0 – любой вектор, в том числе вектор из нулей.
Отсюда получаем при помощи обратной матрицы:
Y*(0) = ∙ ,
Тогда итоговая формула для начала вычислений методом прогонки С.К.Годунова имеет вид:
Y(0) = М∙с + ∙ .
продолжение
--PAGE_BREAK--8 Второй алгоритм для начала счета методом прогонки С.К.Годунова
Этот алгоритм обсчитан на компьютерах в кандидатской диссертации.
Этот алгоритм требует дополнения матрицы краевых условий Uдо квадратной невырожденной:
Начальные значения Y(0), Y(0), Y(0), Y(0), Y*(0) находятся из решения следующих систем линейных алгебраических уравнений:
∙ Y*(0) = ,
∙ Y(0) = , где i = , , , ,
где 0 – вектор из нулей размерности 4х1.
9 Замена метода численного интегрирования Рунге-Кутта в методе прогонки С.К.Годунова
Эта замена формул Рунге-Кутта на формулу теории матриц обсчитана на компьютерах в кандидатской диссертации.
В методе С.К.Годунова как показано выше решение ищется в виде:
Y(x) = Y(x) ∙ c + Y*(x).
На каждом конкретном участке метода прогонки С.К.Годунова между точками ортогонализации можно вместо метода Рунге-Кутта пользоваться теорией матриц и выполнять расчет через матрицу Коши:
Y(x) = K(x— x) ∙Y(x).
Так выполнять вычисления быстрее, особенно для дифференциальных уравнений с постоянными коэффициентами.
И аналогично через теорию матриц можно вычислять и вектор Y*(x) частного решения неоднородной системы дифференциальных уравнений. Или для этого вектора отдельно можно использовать метод Рунге-Кутта, то есть можно комбинировать теорию матриц и метод Рунге-Кутта.
10 Метод половины констант
Этот метод пока не обсчитан на компьютерах.
Выше было показано, что решение системы линейных обыкновенных дифференциальных уравнений можно искать в виде только с половиной возможных векторов и констант. Была приведена формула для начала вычислений:
Y(0) = М∙с + ∙ .
Из теории матриц известно, что если матрица ортонормирована, то её обратная матрица есть её транспонированная матрица. Тогда последняя формула приобретает вид:
Y(0) = М∙с + U∙u
или
Y(0) = U∙u + М∙с
или
Y(0) = ∙ ,
Таким образом записана в матричном виде формула для начала счета с левого края, когда на левом крае удовлетворены краевые условия.
Далее запишем V∙Y(1) = v и Y(1) = K(1←0) ∙Y(0) + Y*(1←0) совместно:
V∙ [ K(1←0) ∙Y(0) + Y*(1←0) ] = v
V∙ K(1←0) ∙Y(0) = v— V∙Y*(1←0)
и подставим в эту формулу выражение для Y(0):
V∙ K(1←0) ∙ ∙ = v— V∙Y*(1←0).
V∙ K(1←0) ∙ ∙ = p.
Таким образом мы получили выражение вида:
D∙ = p,
где матрица D имеет размерность 4х8 и может быть естественно представлена в виде двух квадратных блоков размерности 4х4:
∙ = p.
Тогда можем записать:
D1∙ u + D2 ∙ c= p.
Отсюда получаем, что:
c= D2 ∙ ( p— D1∙ u )
Таким образом, искомые константы найдены.
Далее показано как применять этот метод для решения «жестких» краевых задач.
Запишем
V∙ K(1←0) ∙ ∙ = p.
совместно с K(1←0) = K(1←x2) ∙ K(x2←x1) ∙ K(x1←0) и получим:
V∙ K(1←x2) ∙ K(x2←x1) ∙ K(x1←0) ∙ ∙ = p.
Эту систему линейных алгебраических уравнений можно представить в виде:
[ V∙ K(1←x2) ] ∙ { K(x2←x1) ∙ K(x1←0) ∙ ∙ } = p.
[ матрица ] ∙ { вектор } = вектор
Эту группу линейных алгебраических уравнений можно подвергнуть построчному ортонормированию, которое сделает строчки [матрицы] ортонормированными, {вектор} затронут не будет, а вектор получит преобразование. То есть получим:
[ V∙ K(1←x2) ] {K(x2←x1) ∙ K(x1←0) ∙ } =p.
И так далее.
В итоге поочередного вычленений матриц слева из вектора и ортонормирования получим систему:
D ∙ = p,
Отсюда получаем, что:
c= D2 ∙ (p — D1∙ u)
Таким образом, искомые константы найдены.
11 Применяемые формулы ортонормирования
Эти формулы обсчитаны в кандидатской диссертации.
Взято из: Березин И.С., Жидков Н.П. Методы вычислений, том II, Государственное издательство физико-математической литературы, Москва, 1962 г. 635 стр.
Пусть дана система линейных алгебраических уравнений порядка n:
А=.
Здесь над векторами поставим черточки вместо их обозначения жирным шрифтом.
Будем рассматривать строки матрицы А системы как векторы:
=(,,…,).
Ортонормируем эту систему векторов.
Первое уравнение системы А= делим на .
При этом получим:
++…+=, =(,,…,),
где =, =, =1.
Второе уравнение системы заменяется на:
++…+=, =(,,…,),
где =, =,
=-(,), =-(,).
Аналогично поступаем дальше. Уравнение с номером i примет вид:
++…+=, =(,,…,),
где =, =,
=-(,)-(,)-…-(,),
=-(,)-(,)-…-(,).
Процесс будет осуществим, если система линейных алгебраических уравнений линейно независима.
В результате мы придем к новой системе С=, где матрица С будет с ортонормированными строками, то есть обладает свойством С*С= E, где Е – это единичная матрица.
(Таким образом, решение системы можно записать в виде = С.)
продолжение
--PAGE_BREAK--