bigpo.ru
добавить свой файл
1 2 3







Рассмотрим систему из n нелинейных уравнений, содержащую скалярный параметр q:

(1)

где - достаточно гладкие функции по совокупности аргументов в области их определения. Считая, что вектор с компонентами а вектор с компонентами запишем систему (1) в виде векторного уравнения:

. (2)

Основу численных методов исследования систем нелинейных уравнений составляет теорема о неявной функции. Приведем без доказательства формулировку теоремы применительно к системе (2).

Теорема о неявной функции. Пусть дана система (2), где - вектор-функция векторного аргумента x и скалярного аргумента q, непрерывно дифференцируемая в (n+1)-мерном параллелепипеде с центром в точке , в которой вектор-функция обращается в 0: Если при этом матрица производных не вырождена, то существует окрестность точки где система (2) определяет однозначную непрерывно дифференцируемую вектор-функцию принимающую значение при

Обозначим через D матрицу производных



Если в окрестности решения (2)

, (3)

то тогда по теореме о неявной функции решение x = x(q) системы нелинейных уравнений (2) будет представлено в пространстве [x, q] гладкой пространственной кривой S.

Отметим, что при выполнении условия (3) все аргументы вектор-функции f(x,q), т.е. компоненты векторного аргумента и становятся равноправными.


Перепишем систему (1) к виду

(4)

где

(5)

Введем обозначения:



(6)



Равноправие аргументов проявляется в следующем. Пусть точка принадлежит кривой S, т.е



Если в этой точке выполнено условие

(7)

то по теореме о неявной функции можно считать, что в окрестности система (4) определяет дифференцируемое решение

(8)

т.е. в окрестности играет роль параметра системы (4). Отмечая это обстоятельство, будем записывать (4) в виде, аналогичном (2):

. (9)

Пусть - решение (9), т.е. . Дифференцированием тождества по параметру xk получим систему линейных алгебраических уравнений относительно производных решения (9) по параметру xk:

(10)

Обратим внимание на то, что матрица системы (10) становится известной после завершения итераций по методу Ньютона, который применялся для решения системы (9).

В силу (7),(10) вектор-функцию можно рассматривать также как решение задачи Коши для системы обыкновенных дифференциальных уравнений не разрешенных относительно производных:

. (11)

Используя формулы Крамера, представим (11) в виде:



(12)

Здесь , - определители матриц , в которых -ый столбец заменен на вектор правых частей . Поэтому

(13)

Так как , то можно утверждать, что система (4), где -параметр, определяет S, как интегральную кривую системы дифференциальных уравнений (11), проходящей через точку .

Пусть, в частности,



При этом задача Коши (12), где k = n+1, с учетом (5) принимает вид:



(14)

В тех точках интегральной кривой (14), в которых при q = q*, графики функций xi = xi(q) имеют вертикальную касательную. Таким образом, по отношению к (14) условие определяет особую точку типа «поворот», в которой происходит бифуркация числа (изменение числа решений).

Возможно, что в некоторых точках кривой S ранг матрицы D окажется меньше n. Такие точки называются существенно особыми. Предлагаемый метод продолжения решения по параметру позволяет численно построить дугу гладкой пространственной кривой S между двумя особыми точками, если особые точки имеют место.

Приведем описание одного регулярного шага построения кривой S в окрестности точки методом продолжением по параметру решения системы (4) с использованием так называемой «параметризации». При этом предполагается, что роль параметра в (4) играет компонента xk, а система линейных алгебраических уравнений



определяет в точке вектор производных .

Параметризацией называется правило выбора компоненты μ = xj, в точке кривой S, которая в качестве параметра системы (4) используется для продолжения решения (4) на один шаг в окрестности . При этом производные решения (4) по параметру μ в точке определяются нормировкой компонент вектора производных . Так как параметр μ выбирается на один шаг, то он называется «текущим параметром».

В связи с выбором текущего параметра μ компоненты вектора производных подвергаются нормировке. Для этого определяется максимальная по модулю компонента. Пусть



Тогда, как следует из (12) с учетом (13), имеем



т.е.


        1. Так как при выполнении условия (15) то в окрестности точки компонента xj также может быть параметром системы (4). В этом случае система (4) определяет дифференцируемое решение




Если выполнено условие (15), то производные решения по параметру xj в точке находятся нормировкой производных по xk:



При этом модули компонент (16) вектора производных по xj ограничены единицей. Если условие (15) не выполняется, то нормировка не производится, так как в этом случае модули компонент вектора производных и так ограничены единицей. Параметром системы (4) остается xk. Таким образом, выбор параметра μ в точке происходит по правилу:

(17)

где ω определено в (10). Для производных решения по параметру μ соответственно имеем:





(Здесь содержится тривиальное выражение: .)

Смысл проделанных преобразований состоит в том, что, зная решение и производную решения по текущему параметру, можно использовать их для задания начального приближения решения системы (4) в достаточно малой окрестности точки σ0. Кроме того, по непрерывности модули компонент вектора производных по текущему параметру будут мало отличаться от 1 в достаточно малой окрестности точки σ0 .

Воспользуемся обозначениями, аналогичными (6). Пусть m – индекс текущего параметра μ, определенного при параметризации в точке σ0 согласно (17), где μ = μ0. При этом



который является (по аналогии с (9)), решение системы

.

Если μ = μ0, то в этих обозначениях вектор Xμ с компонентами будет решением системы



Обозначим, далее, через шаг по текущему параметру μ. Тогда, например, формулы

(18)

могут быть использованы в методе Ньютона для задания начального приближения решения системы:

. (19)

Величина шага hμ по параметру μ в формулах (18),(19) может варьироваться в зависимости от результатов сходимости итераций по Ньютону при решении системы (19). Первоначально hμ присваивается значение пробного шага и по формулам (18) вычисляется начальное приближение решения (19). Если при этом число итераций достигает ограничения, но норма невязок на итерациях все еще больше заданного числа, то шаг hμ делится пополам и снова задается начальное приближение (18). Вычисления прекращаются, если величина шага становится меньше допустимого. В противном случае считается, что итерационный процесс сошелся, т.е. получено решение (19) и, следовательно, становятся известными координаты точки , где , соседней с точкой σ0. Продолжение решения (4) на один шаг из точки σ0 в точку σ1 завершено.

В точке σ1 выполняются те же действия, что и в точке σ0. А именно.Применением параметризации определяется текущий параметр λ, значение которого при равно . Это означает, что найдено решение системы



а также вектор производных решения по параметру λ в точке σ1, модули компонент которого ограничены единицей. В результате появляется возможность задания начального приближения решения системы



где hλ – пробный шаг по текущему параметру λ, и т.д.

Замечание. Если точка σ0 – начало кривой S (стартовая позиция), то начальное приближение берется в виде (18), где пробный шаг задается.


Если точка σ0 отвечает текущей позиции, то пробный шаг hλ в точке σ1 при задании начального приближения решения системы



берется равным , где с > 1 – «параметр разгона».

Отметим, что в текущей позиции начальное приближение может задаваться в виде кубической параболы с привлечением производных решения по λ при и . Для вычисления производных при достаточно разделить производные решения по параметру при на производную . При этом начальное приближение имеет вид:





где

,

.

Процесс построения кривой S завершается на границе задаваемой области изменения параметра q. В результате имеем таблицу значений используя которую можно построить зависимость компонент вектора Xk от xk , где клюбой столбец этой таблицы, в том числе и (n+1), содержащий значения параметра q. График зависимости x = x(q) является итогом численного исследования системы (2).

Иллюстрационный пример построения гладкой пространственной кривой S для системы



.

представлен на рис.1, 2. Графики проекций S указывают на существование множественности решений рассматриваемой системы уравнений.




РИС.1.Пример гладкой пространственной кривой S, которая представляет решение системы из 2-х уравнений с параметром q:

.




РИС.2 Диаграмма множественности решений в виде графиков компонент и .

Применение метода продолжения решения по параметру носит характер численного эксперимента, в котором присутствует ряд вычислительных проблем, требующих индивидуального подхода в зависимости от конкретного вида каждой из рассматриваемых нелинейных систем уравнений (как и при численном исследовании нелинейных проблем вообще). Отметим некоторые из них.

  1. Задание начального приближения в методе Ньютона в стартовой позиции.

  2. Определение области изменения параметра системы. Вообще говоря, если известно решение системы при некотором значении параметра , то можно лишь утверждать, что по непрерывности решение существует лишь в некоторой окрестности , размеры которой заранее неизвестны и устанавливаются в ходе численного эксперимента.

  3. Подбор параметров метода продолжения решения (параметры метода Ньютона, параметры адаптации текущих шагов), от которых может зависеть вид численно построенной пространственной кривой S . Отметим, что в процессе продолжения находится последовательность точек, принадлежащих графику. Поэтому для правильного представления графика, необходимо согласование параметров метода с гладкостью S.

Метод продолжения решения по параметру является инструментом численного исследования, результаты которого тем достовернее, чем полнее априорные представления о свойствах изучаемой системы нелинейной системы.

      1. Тема 2. Продолжение решения по параметру, как задача Коши.



Остановимся на варианте продолжения решения системы нелинейных уравнений



в котором проблема формулируется как задача Коши (14):



Здесь





- k-я точка интегральной кривой.

Требуется построить интегральную кривую, проходящую через точку . Предполагается, что известны формулы метода численного интегрирования на один шаг и что для интегрирования выбран метод с автоматическим выбором шага.

Пусть - текущий параметр в точке , где . При этом для выбора текущего параметра применялась параметризация. В результате использования формул интегрирования на один шаг по

задачи Коши



становится известной точка , где применением параметризации выбирается текущий параметр . И так далее.

Пусть i = 1,2,…, n+1, где xn+1 = q, - известное решение (1) при некотором значении параметра q, а - текущий параметр. Тогда можно считать, что это решение является результатом интегрирования на один шаг системы обыкновенных дифференциальных уравнений с независимым аргументом :



где вектор с компонентами ; матрица производных со столбцами . Поскольку производные решения по параметру известны, то это дает возможность применением параметризации определить новый текущий параметр xj =. При этом построение интегральной кривой на один шаг осуществляется интегрированием системы дифференциальных уравнений



где вектор с компонентами ; матрица производных со столбцами . И т.д.

      1. Тема 3. Метод Кубичека.



В методе Кубичека параметром для продолжения решения системы нелинейных уравнений (2) является длина дуги s гладкой пространственной кривой S. Очевидно, в параметрическом представлении решения (2) в виде имеет место однозначная зависимость x и q от s.

Пусть - произвольная точка S,

(20)

Дифференцированием (20) по s получим:



или

(21)

где nх(n+1) – матрица, ранг которой равен n. Дополним систему уравнений (12) выражением элемента длины дуги:

(22)

В стартовой точке продолжения выполняются условия:

при (23)

Воспользуемся методом Гаусса с выбором главного элемента по всей матрице для преобразования (21) к виду:

(24)

где B – верхняя треугольная матрица, b – преобразованный k-ый столбец матрицы , компоненты вектора образуются из компонент вектора в результате перестановок, связанных с методом Гаусса.


Из (24) следует система линейных алгебраических уравнений с невырожденной матрицей B относительно производных компонент вектора u по параметру μ

(25)

Поскольку , то играет роль текущего параметра. При этом элемент дуги записывается в виде:

(26)

Выбор знака согласуется с условием возрастания длины дуги.

Из тождества

,

а также уравнений (25),(26) следуют равенства, которые можно рассматривать как систему дифференциальных уравнений относительно u = u(s) и μ = μ(s). Имеем:

, . (27)

Для вычисления правых частей (27) необходимо найти решение системы линейных алгебраических уравнений (28) с верхней треугольной матрицей B. Напомним, что соответствие между компонентами векторов (x(s), q(s)) и (u(s), μ(s)) устанавливается перестановками метода Гаусса. Начальные условия следуют из (28). Тем самым задача продолжения решения системы (2) по длине дуги пространственной кривой S сведена к задаче Коши для системы (27) с независимым аргументом s .

Приведем схему вычислений по методу Кубичека.

1. Предполагается, что известен способ решения системы нелинейных уравнений (2) при q = q0:



2. При q = q0 вычисляются элементы матрицы производных

.

3. Матрица G применением метода Гаусса с выбором главного элемента по всей матрице преобразуется к матрице Q = [B, b].

В ходе преобразования запоминаются перестановки столбцов матрицы G. При этом столбец с номером k матрицы G занимает положение столбца b матрицы Q, а μ = xk становится текущим параметром. Элементы вектора (x, q) после перестановок становятся элементами вектора (u, μ).

4. Решается система линейных уравнений



с верхней треугольной матрицей B относительно z = du/.

5. Вычисляются правые части системы дифференциальных уравнений (27):

, .

Знак правых частей согласуется с условием возрастания длины дуги,

,

с учетом знака приращения текущего параметра.

6. Выполняется шаг интегрирования, при котором находится решение (27) (u1, μ1) по формулам выбранного метода интегрирования. С использованием известных перестановок восстанавливаются элементы вектора (x1, q1)

7. При q = q1 вычисляются элементы матрицы производных



8. Выполняется пункт 3 и т.д.

Сделаем ряд замечаний к методу Кубичека. В методах продолжения решения системы нелинейных уравнений по параметру, сформулированных как задача Коши, указывается способ вычисления правых частей системы дифференциальных уравнений, что формально позволяет использовать одну из стандартных программ интегрирования для численного построения решения x = x(q) системы нелинейных уравнений (2). Однако, возможны случаи, когда система дифференциальных уравнений оказывается жесткой. Возникающие в связи с этим вычислительные проблемы (например, применение явного метода интегрирования из-за сложности вычисления производных правых частей) обусловлены только выбором метода решения (2).
      1. Раздел 4.

      2. ЧИСЛЕННЫОЕ ИССЛЕДОВАНИЕ НЕЛИНЕЙНЫХ КРАЕВЫХ ЗАДАЧ. МЕТОД ПРОДОЛЖЕНИЯ ПО ПАРАМЕТРУ.




      1. Тема 1. Продолжение решения по параметру в методе множественной стрельбы.



Рассмотрим нелинейную краевую lдля системы из n дифференциальных уравнений:

(1)

где y, f и g – векторы с компонентами yk, fk и gk, k = 1,2,…,n, соответственно, q – скалярный параметр. Вектор-функции f(x,y,q) и g(y(a),y(b),q) предполагаются достаточно гладкими по совокупности их аргументов. Предполагается, что краевая задача (1) хорошо обусловлена и имеет решение y = y(x,q), дифференцируемое как по x, так и по

По определению также хорошо обусловленной будет краевая задача относительно производной решения (1) по параметру q, которую обозначим как вектор-функцию v(x):

.

Подставим в (1) решение y = y(x,q). После дифференцирования полученных тождеств по q приходим к выводу, что v(x) является решением линейной краевой задачи:

(2)

где

,

(3)

EMBED Equation.3

Для удобства здесь использованы обозначения:

Совместное решение краевых задач (1) и (2) позволяет реализовать продолжение решения (1) на основе метода множественной стрельбы.

Напомним (раздел 2, тема 3), что метод множественной стрельбы применительно к нелинейной краевой задаче

(4)

сводится к заданию сетки,

(5)

и решению методом Ньютона системы нелинейных уравнений относительно сеточных значений искомой вектор-функции y(x). Вектор-функция и матрица производных определяются в ходе итераций на решениях известной серии задач Коши.

Очевидно, в случае (1) аналогичная система нелинейных уравнений относительно сеточных значений решения y = y(x,q) будет иметь вид:

. (6)

Предполагается, что в некоторой области изменения параметра q система (4) определяет зависимость P = P(q), которая в пространстве (P, q) имеет представление в виде гладкой пространственной кривой Σ.

Пусть det(ΨP) ≠ 0 в некоторой точке кривой Σ. Тогда производная решения (6) по параметру находится из системы линейных алгебраических уравнений:

(7)

Отметим, что матрица ΨP вычисляется на итерациях по Ньютону при решении системы (6).

Как известно (раздел 3), применение к системе (6) метода продолжения по параметру для построения кривой Σ , а вместе с этим и решения краевой задачи (1) в зависимости от параметра, требует вычисления матрицы производных

(8)

Приведем формулировку серии задач Коши, связанную с вычислением матрицы Ω.

Как и в методе множественной стрельбы для решения краевой задачи (4), сформулируем серию задач Коши, которая с учетом краевых условий (1) определяет вектор-функцию и матричную функцию . Имеем:



(9)



где pk, k = 1,2,…,m+1, - векторная компонента вектора P.

Пусть y=y(x,p[k],q) и Y=Y(x,p[k],q) – решения (9) на отрезке [xk,xk+1]. При этом из условий непрерывности решения краевой задачи (1) и краевых условий следуют формулы для вычисления векторных компонент ψ[j] , j =0, 1,…,m, вектора Ψ:



. (10)

Матрица производных имеет вид:



Определим теперь вектор-столбец Ψq матрицы Ω. C этой целью воспользуемся методом множественной стрельбы для численного решения линейной краевой задачи (2). Введем обозначения для сеточных значений вектор-функции v(x):



При этом вектор является векторной компонентой составного вектора dP/dq, определяемого из системы линейных алгебраических уравнений (7). В соответствии с методом множественной стрельбы (раздел 2, тема 1) рассматривается серия задач Коши:




(12)




следующая страница >>