Учебное пособие: Методы решения краевых задач, в том числе "жестких" краевых задач

Как выполнять построчное ортонормирование систем линейных алгебраических уравнений можно посмотреть в [Березин, Жидков].

Дополним прямоугольную горизонтальную матрицу 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 + Yx + Yx/2!  +  …,  где Y=Y(0), Y= Y(0), …       (2).

Из (1) почленным дифференцированием при А=const и F(x)=0 получим:

Y= AY= AY,        Y= A Y = AY,     (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)= eYa(x).      (7)

Продиффренцируем (7) и подставим в (1). Получим:

eYa(x) = F(x).       (8)

При получении (8) учитывалось, что:

 =  = A + A x + A x/2!  + … = A e.

Из (8) следует, что:

Ya(x) = c + .      (9)

Подставим в (7) и получаем:

Y(x) = ec + 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. Метод для численного интегрирования дифференциальных уравнений.

Читали нам как-то в бауманке численные методы решения дифференциальных уравнений. И, кажется, приводили аналитический вывод формул одного из авторов. Или это просто мелькнуло в учебнике (я имею в виду вывод формул). Уже не очень помню. Запомнилась только собственная мысль, что людям вообще-то проще всего даются геометрические аналогии и выводы, сделанные на основе понятных геометрических картинок. Ну, вот тогда я и нарисовал один из вариантов численного решения дифференциальных уравнений и помню даже перевёл геометрические картинки в буквенные формулы приближённых вычислений. Сейчас повторно выводить буквенные формулы для численного интегрирования дифференциальных уравнений мне не кажется интересным. А вот привести картинки тех студенческих мыслей вполне можно для обсуждения.

Далее идёт картинка с текстом и с рисунками численного интегрирования: