Методы решения краевых задач, в том числе «жестких» краевыхзадач
Методы Алексея Юрьевича Виноградова
1Введение
Напримере системы дифференциальных уравнений цилиндрической оболочки ракеты –системы обыкновенных дифференциальных уравнений 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 это вектор частного решения неоднороднойсистемы дифференциальных уравнений.
2Случай переменных коэффициентов
Этотвариант рассмотрения переменных коэффициентов проверялся в кандидатскойдиссертации.
Изтеории матриц [Гантмахер] известно свойство перемножаемости матричных экспонент(матриц Коши):
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/>.
3Формула для вычисления вектора частного решения неоднородной системыдифференциальных уравнений
Этаочень простая формула еще не обсчитана на компьютерах. Вместо неё обсчитываласьзначительно ранее выведенная и гораздо более сложная формула, приведенная в:
Численныйметод переноса краевых условий для жестких дифференциальных уравненийстроительной механики Журнал «ММ», Том: 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.
Дляслучая переменных коэффициентов A=A(x) можно использовать прием разделенияучастка (x/> — x/>) интервала интегрирования на малыеподучастки, где на подучастках коэффициенты можно считать постоянными A(x/>)=const и тогда вектор частного решения неоднородной системыдифференциальных уравнений Y*(x/>←x/>) будет на участке складываться из соответствующихвекторов подучастков, на которых матрицы Коши приближенно вычисляются припомощи формул с постоянными матрицами в экспонентах.
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) = />.
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) =М/>∙с + />∙ />.
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, где Е – это единичнаяматрица.
(Такимобразом, решение системы можно записать в виде />=С/>/>.)
12Вывод формул, позаимствованный из «Теории матриц» Гантмахера
Системалинейных обыкновенных дифференциальных уравнений с постоянными коэффициентамиимеет вид:
Y/>(x) = A Y(x) + F(x). (1)
РазложимY(x) в ряд Маклорена по степеням x:
Y(x)=Y/> + Y/>x + Y/>x/>/2! + …, где Y/>=Y(0), Y/>= Y/>(0), … (2).
Из(1) почленным дифференцированием при А=const и F(x)=0 получим:
Y/>= AY/>= A/>Y, Y/>= A Y/> = A/>Y, (3)
Положивв (3) x=0 и подставив в (2) получим:
Y(x) = Y/> + Ax Y/> + A/> x/>/2! Y/> + … = e/> Y/>, (4)
гдеe/> = E + Ax + A/> x/>/2! + …, где Е – единичная матрица.(5)
Еслипринять x=x/>, то (4) заменится на
Y(x) = e/> Y(x/>), (6)
Рассмотримслучай A=const и F≠0.
Введемв рассмотрение вектор-функцию Ya(x) в виде: Y(x)= e/>Ya(x). (7)
Продиффренцируем(7) и подставим в (1). Получим:
e/>Ya/>(x) = F(x). (8)
Приполучении (8) учитывалось, что:
/> = /> = A + A/> x + A/> x/>/2! + … = A e/>.
Из(8) следует, что:
Ya(x) = c + />. (9)
Подставимв (7) и получаем:
Y(x)= e/>c + e/>/>. (10)
Положивx=x/> в (10) получим:
c = e/> Y(x/>). (11)
Окончательнополучаем:
Y(x)= e/> Y(x/>) + e/>/>. (12)
Мойотец предложил использовать и другую (гораздо более эффективную по временисчета) матричную формулу вместо матричной экспоненты – что-то на основе Вольтерра.Это есть в статье в журнале «Математическое моделирование»:
Численныйметод переноса краевых условий для жестких дифференциальных уравненийстроительной механики Журнал «ММ», Том: 14 (2002), Номер: 9, 3 стр. 1409-003r.pdf
Ятут как-то на днях поговорил с отцом и вот что любопытное насчет методавыяснилось.
Когдая сам решал своим методом «переноса краевых условий» «жесткие» краевые задачи,то я в каждой рассматриваемой точке x=x* решал соответствующую систему линейныхалгебраических уравнений для нахождения решения Y(x*) в каждой рассматриваемойточке x=x*.
Амой отец утверждает, что «жесткими» бывают только краевые задачи, а начальныезадачи «жесткими» не бывают. И поэтому он находил решение Y(x*) в какой-тоодной точке x=x* по моему методу, а дальше решал (влево и вправо отрассматриваемой точки x=x*) как задачу Коши (от найденных начальных условийY(x*) в этой одной точке): Y(x)=K(х←x*)·Y(x*)+Y*(x←x*). И если это так, то так решать,естественно, гораздо быстрее, так как надо только домножать на матрицу Коши(матрициант, матричную экспоненту) вместо решения систем линейныхалгебраических уравнений.
Ядумаю, что это очень существенное уточнение по скорости счета.
P.P.P.S. Метод решения«жестких» начальных задач, то есть «жестких» задач Коши. Придуман вечером 16января 2008 года.
Методпереноса краевых условий работает успешно. Следовательно, по аналогии можнопытаться решать «жесткие» начальные задачи, то есть «жесткие» задачи Коши.
Пустьдана начальная задача:
Y/>(x) = A·Y(x), Y(0) = Yнач,
Можемзаписать:
Y(0)= K(0←x) · Y(x),
K(0←x)· Y(x) = Yнач
K(0←x3)· K(x3←x2) · K(x2←x) · Y(x) = Yнач
[ K(0←x3) ] · { K(x3←x2) ·K(x2←x) · Y(x) } = Yнач
[матрица ] { вектор } вектор
Выполняемпострочное ортонормирование этой системы линейных алгебраических уравнений сквадратной матрицей коэффициентов и получаем систему с ортонормированнымистрочками в квадратной матрице:
[K(0←x3) ]орт · {K(x3←x2) · K(x2←x) · Y(x)} = Yнач_орт
Аналогичнозаписываем
[[K(0←x3) ]орт · K(x3←x2) ] ·{ K(x2←x) · Y(x) } = Yнач
[ матрица ] { вектор } вектор
Далеевыполняем построчное ортонормирование и получаем:
[[K(0←x3) ]орт · K(x3←x2)]орт · {K(x2←x) · Y(x)} = Yнач_2орт
Аналогичнополучаем
[[K(0←x3) ]орт · K(x3←x2)]орт · K(x2←x)]орт · Y(x) = Yнач_3орт
Тоесть получили итоговую систему линейных алгебраических уравнений для нахождениявектора Y(x). Для этого вектор Yнач_3орт надо домножить слева на матрицутранспонированную по отношению к левой ортонормированной матрице, так как еслиматрица ортонормированна, то обратная ей есть транспонированная матрица.
P.P.P.P.S. 11 сентября 2009:
Долгобыло лень записывать, но вдруг кто-то сам не сразу догадается – ещё один методрешения «жестких» начальных задач, то есть «жестких» задач Коши. Хотя у меняесть такое подозрение, что «жесткими» бывают только краевые задачи, а начальныезадачи «жесткими» НЕ бывают. Но на всякий случай приведу метод решения«жестких» начальных задач.
Начальныеусловия имеют вид:
Y(0) =Yнач.
Полноерешение системы обыкновенных дифференциальных уравнений Y/>(x) = A · Y(x) + F(x) имеет вид:
Y(x) = K(x←0) · Y(0) + Y*(x←0).
Илиможно записать:
Y(0) =K(0←x1) · Y(x1)+ Y*(0←x1).
Подставимэто выражение в краевые условия и получим:
K(0←x1) · Y(x1) + Y*(0←x1) = Yнач
или
K(0←x1) · Y(x1) = Yнач — Y*(0←x1)
или
K1 · Y(x1) = Y1.
Проортонормируемэто выражение построчно и получим эквивалентное выражение:
K1орто· Y(x1) = Y1орто.
Тогда
Y(x1) = (K1орто)транспонир · Y1орто.
Подставимвместо Y(x1) выражение через Y(x2) и получим:
K(x1←x2) · Y(x2)+ Y*(x1←x2) =(K1орто)транспонир · Y1орто
или
K(x1←x2) · Y(x2)= (K1орто)транспонир · Y1орто — Y*(x1←x2)
или
K2 · Y(x2) = Y2.
Проортонормируемпострочно и получим эквивалентное выражение:
K2орто· Y(x2) = Y2орто.
Тогда:
Y(x2) = (K2орто)транспонир · Y2орто.
Итак далее.
P.P.P.P.P.S. Метод для численного интегрирования дифференциальныхуравнений.
Читалинам как-то в бауманке численные методы решения дифференциальных уравнений. И,кажется, приводили аналитический вывод формул одного из авторов. Или это простомелькнуло в учебнике (я имею в виду вывод формул). Уже не очень помню.Запомнилась только собственная мысль, что людям вообще-то проще всего даютсягеометрические аналогии и выводы, сделанные на основе понятных геометрическихкартинок. Ну, вот тогда я и нарисовал один из вариантов численного решениядифференциальных уравнений и помню даже перевёл геометрические картинки вбуквенные формулы приближённых вычислений. Сейчас повторно выводить буквенныеформулы для численного интегрирования дифференциальных уравнений мне не кажетсяинтересным. А вот привести картинки тех студенческих мыслей вполне можно дляобсуждения.
Далееидёт картинка с текстом и с рисунками численного интегрирования:
/>