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

• СОДЕРЖАНИЕ
• 1. Методы решения систем нелинейных уравнений
• 2. Векторная запись нелинейных систем
• 3. Метод Ньютона, его реализации и модификации
• 3.1 Метод Ньютона
• 3.2 Модифицированный метод Ньютона
• 3.3 Метод Ньютона с последовательной аппроксимацией матриц
• 3.4 Разностный метод Ньютона
• 3.5 Обобщение полюсного метода Ньютона на многомерный случай
• 4. Численный пример
• 5. Реализация метода Ньютона в среде MATLAB
• Список используемой литературы

1. Методы решения систем нелинейных уравнений

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

2. Векторная запись нелинейных систем.

Пусть требуется решить систему уравнений

(2.1)

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

,

данную систему (2.1) можно записать одним уравнением

(2.1а)

относительно векторной функции F векторного аргумента х. Таким образом, исходную задачу можно рассматривать как задачу о нулях нелинейного отображения. В этой постановке она является прямым обобщением основной задачи — построения методов нахождения нулей одномерных нелинейных отображений. Фактически это та же задача, только в пространствах большей размерности. В любом случае следует позаботиться о правомочности тех или иных операций над векторными переменными и векторными функциями, а также о сходимости получаемых таким способом итерационных процессов. Часто теоремы сходимости для этих процессов являются тривиальными обобщениями соответствующих результатов, полученных для методов решения скалярных уравнений. Однако не все результаты и не все методы можно перенести со случая п = 1 на случай п ?2. Например, здесь уже не будут работать методы дихотомии, поскольку множество векторов не упорядочено. В то же время, переход от n = 1 до n?2 вносит в задачу нахождения нулей нелинейного отображения свою специфику, учет которой приводит к новым методам и к различным модификациям уже имеющихся. В частности, большая вариативность методов решения нелинейных систем связана с разнообразием способов, которыми можно решать линейные алгебраические задачи, возникающие при пошаговой линеаризации данной нелинейной вектор-функции F (x).

3. Метод Ньютона, его реализации и модификации

3.1 Метод Ньютона

Пусть () — некоторая последовательность невырожденных вещественных n x n-матриц. Тогда, очевидно, последовательность задач

k = 0,1,2,…

имеет те же решения, что и исходное уравнение (2.1а), и для приближенного нахождения этих решений можно формально записать итерационный процесс

k = 0,1,2,… (3.1.1)

при. В случае — это действительно МПИ с линейной сходимостью последовательности () Если же различны при разных k, то формула (3.1.1) определяет большое семейство итерационных методов с матричными параметрами. Рассмотрим некоторые из методов этого семейства.

Положим, где

— матрица Якоби вектор-функции F (x). Подставив это в (3.1.1), получаем явную формулу метода Ньютона

(3.1.2)

обобщающего на многомерный случай скалярный метод Ньютона (5.14). Эту формулу, требующую обращения матриц на каждой итерации, можно переписать в неявном виде:

. (3.1.3)

Применение (3.1.3) предполагает при каждом k = 0,1,2,… решение линейной алгебраической системы

относительно векторной поправки, а затем прибавление этой поправки к текущему приближению для получения следующего:

.

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

Сравнивая (3.1.3) с формальным разложением F (x) в ряд Тейлора

видим, что последовательность () в методе Ньютона получается в результате подмены при каждом k=0,1,2,… нелинейного уравнения F (x) = 0 или, что-то же (при достаточной гладкости F (x)), уравнения

линейным уравнением

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

Новым, по сравнению со скалярным случаем, фактором, осложняющим применение метода Ньютона к решению n-мерных систем, является необходимость решения n-мерных линейных задач на каждой итерации (обращения матриц в (3.1.2) или решения СЛАУ в (3.1.3)), вычислительные затраты на которые растут с ростом n, вообще говоря, непропорционально быстро. Уменьшение таких затрат — одно из направлений модификации метода Ньютона.

3.2 Модифицированный метод Ньютона

Если матрицу Якоби F'(х) вычислить и обратить лишь один раз — в начальной точке, то от метода Ньютона (3.1.2) придем к модифицированному методу Ньютона

(3.2.1)

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

Компромиссный вариант — это вычисление и обращение матриц Якоби не на каждом итерационном шаге, а через несколько шагов (иногда такие методы называют рекурсивными).

Например, простое чередование основного (3.1.2) и модифицированного (3.2.1) методов Ньютона приводит к итерационной формуле

(3.2.2)

где k = 0,1,2,… За здесь принимается результат последовательного применения одного шага основного, а затем одного шага модифицированного метода, т. е. двухступенчатого процесса

(3.2.3)

Доказано, что такой процесс при определенных условиях порождает кубически сходящуюся последовательность ().

3.3 Метод Ньютона с последовательной аппроксимацией матриц

Задачу обращения матриц Якоби на каждом k-м шаге метода Ньютона (3.1.2) можно попытаться решать не точно, а приближенно. Для этого можно применить, например, итерационный процесс Шульца, ограничиваясь минимумом — всего одним шагом процесса второго порядка, в котором за начальную матрицу принимается матрица, полученная в результате предыдущего (k-1)-го шага. Таким образом приходим к методу Ньютона с последовательной аппроксимацией обратных матриц:

(3.3.1)

где k = 0,1,2,…, а и — начальные вектор и матрица (). Этот метод (будем называть его более коротко ААМН — аппроксимаиионный аналог метода Ньютона) имеет простую схему вычислений — поочередное выполнение векторных в первой строке и матричных во второй строке его записи (3.3.1) операций. Скорость его сходимости почти так же высока, как и у метода Ньютона. Последовательность () может квадратично сходиться к решению уравнения F (x)=0 (при этом матричная последовательность () также квадратично сходится к, т. е. в нормально развивающемся итерационном процессе (3.3.1) должна наблюдаться достаточно быстрая сходимость () к нулю).

Применение той же последовательной аппроксимации обратных матриц к простейшему рекурсивному методу Ньютона (3.2.2) или, что-то же, к двухступенчатому процессу (3.2.3) определяет его аппроксимацирнный аналог

(3.3.2)

как и (3.2.2), также можно отнести кметодам третьего порядка. Доказательство кубической сходимости этого метода требует уже более жестких ограничений на свойства F (х) и близость к, к, чем в предыдущем методе. Заметим, что к улучшению сходимости здесь может привести повышение порядка аппроксимации обратных матриц, например, за счет добавления еще одного слагаемого в формуле для подсчета :

3.4 Разностный метод Ньютона

На базе метода Ньютона (3.1.2) можно построить близкий к нему по поведению итерационный процесс, не требующий вычисления производных. Сделаем это, заменив частные производные в матрице Якоби J (x) разностными отношениями, т. е. подставив в формулу (3.1.1) вместо матрицу где

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

3.5 Обобщение полюсного метода Ньютона на многомерный случай

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

(3.5.1)

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

(3.5.2)

(вторая компонента — ордината точки пересечения указанных прямых — после вычисления значения может дать информацию об отклонении от функции в точке ее локальной аффинной модели, каковой является проведенная в точке касательная). Ясно, что получаемое из системы (3.5.2) значение тождественно его.

Рассмотренный векторный подход к построению одномерного полюсного метода Ньютона служит ключом для его распространения на двумерный случай на основе таких же геометрических, но уже пространственных соображений.

Пусть требуется найти приближенное решение двумерной нелинейной системы в предположении непрерывной дифференцируемости входящих в нее функций f (x, у) и g (x, у) в некоторой области G, содержащей искомое решение х* =(х*; у*) и приближения к нему k = 0,1,2,…

Будем считать, что уже найдено k-е приближение к решению х* и нужно получить правило перехода к (k+1)-му приближению. В сделанном предположении о гладкости функций f (x, у) и g (x, у) можно провести касательные плоскости в точке () определяемым ими поверхностям

z = f (x, y) и z = g (x, y). (3.5.3)

Эти плоскости задаются текущими векторами

и нормалями

соответственно, т. е. аналогично первому из условий (3.5.1) должно быть

иначе,

. (3.5.4)

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

Чтобы осуществить поставленную цель, зафиксируем в R3 две несовпадающие между собой и с точки — полюсы и. Через указанные три точки можно провести единственную плоскость (которая здесь играет роль прямой, проходящей через полюс и точку (хк; 0) в одномерной ситуации). Взяв текущую точку М (х; у; z) и образовав текущий вектор этой третьей плоскости, можно задать ее условием компланарности трех векторови (что служит аналогом второго из условий (3.5.1)).

Запишем совокупность всех трех описанных средствами векторной алгебры плоскостей в координатной форме. Имеем:

Первые две координаты вектора (x;y;z), служащего решением полученной системы уравнений, считаем искомым приближением ().Введя поправки

(3.5.5)

эту систему превращаем в систему уравнений относительно неизвестных и z:

(3.5.6)

Для исключения вспомогательной переменной z из линейной системы (3.5.6) выразим ее из третьего уравнения. Обозначив

, (3.5.7)

раскрываем фигурирующий в (3.5.6) определитель по элементам первой строки:

Отсюда находим выражение

(3.5.8)

подставляя которое в первые два уравнения системы (3.5.6), приходим к двумерной линейной системе

(3.5.9)

Фактически эта система вместе с обозначениями (3.5.7) и определяет двумерный полюсный метод Ньютона для нелинейной системы. Надя их нее поправки, в соответствии с равенствами (3.5.5) получаем очередное приближение :

.

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

Пусть задана нелинейная система, функции (образующими вектор) в точке, можно описать матрично-векторным уравнением

(3.5.10)

где — n-мерный вектор, каждой компонентой которого служит вспомогательная переменная, входящая в уравнения гиперповерхностей .

Зададим n полюсов (i=1,2,…, n) так, чтобы они не принадлежали одной гиперплоскости пространства. Через все эти полюсы и точку (), определяемую известным приближением к решению системы, проводим гиперплоскость, уравнение которой аналогично двумерному случаю задаем условием равенство нулю определителя (n+1) порядка:

. (3.5.11)

Векторно-матричное уравнение (3.5.10) и скалярное уравнение (3.5.11), в принципе, уже определяют векторный n-полюсный метод Ньютона для построения приближенной к решению системы. Чтобы записать соответствующую линейную систему относительно поправок

(3.5.12)

(аналогичную схеме (3.5.9) двумерного случая), введем следующие обозначения. Положим

,

и образуем квадратную (n+1)-мерную матрицу следующей структуры:

.

Тогда на основе (3.5.10), (3.5.11) имеем (n+1)-мерную систему уравнений относительно неизвестных :

(3.5.13)

Как и в двумерном случае, из второго уравнения этой системы выражаем вспомогательную неизвестную :

(3.5.14)

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

Заменив в (3.5.13) все компоненты вектора z найденным их значением (3.5.14), приходим к следующему линейному векторно матричному уравнению относительно вектора-поправки :

(3.5.15)

Где

. (3.5.16)

Уравнение (3.5.15) вместе со связью (3.5.12), согласно которой

(3.5.17)

является неявной формой пполюсного метода Ньютона для уравнения (2.1а).

Совокупности формул (3.5.15)-(3.5.17) можно придать другой вид:

(3.5.18)

который удобно трактовать как явный метод Ньютона со своеобразной коррекцией матриц Якоби путем прибавления к ним формирующихся по заданному правилу матриц. Как и в одномерном случае, для ускорения сходимости последовательности приближений полюсы целесообразно изменять в такт с изменением значений функций, и в самом простом случае есть смысл фиксировать матрицу С, а вектор брать равным или ;

4. Численный пример

Начальное приближение:

Вектор-функция:

Матрица Якоби вектор-функции:

Вычисляем корень по формуле метода Ньютона c точностью :

Ответ:

нелинейное уравнение решение ньютон

5. Реализация метода Ньютона в среде MATLAB

М-функция.

function nwt = newton (x, y, e, F0, F1,dF0x, dF0y, dF1x, dF1y)

for i = 1:1 000 000

F=[F0(x, y); F1(x, y)];

dF=[dF0x (x, y) dF0y (x, y); dF1x (x, y) dF1y (x, y)];

Z = [x;y] - dF^(-1)*F;

b = sqrt ((x-Z (1))^2+(y-Z (2))^2);

if b > e

x = Z (1);

y = Z (2);

else break;

end

end

disp ('Количество итераций'); disp (i);

nwt = Z;

Функция в качестве входных параметров принимает начальное приближение (), функцию (), частные производные функции и точность e.

Пятая строка находит новую точку приближения. Шестая строка вычисляет норму между текущим и следующим приближением. Строки восемь и девять запоминают точку начального приближения.

Процесс нужно продолжать до тех пор пока. Если, процесс завершить и получим решение .

Теперь напишем скрипт который покажет работу нашей М-функции.

Найдем решение заданной системы нелинейных уравнений

при начальном приближении x=0, y=-1, с точностью до 0.001:

Скрипт.

% Решить систему уравненийу методом Ньютона

% sin (x+y)-1.6x

% x²+y²−1

% Введём функцию F (x)(систему функций)

F0 = inline ('sin (x+y)-1.6*x');

F1 = inline ('x²+y²−1');

disp (F0); disp (F1);

% Их производны

dF0x = inline ('cos (x+y)-1.6');

dF0y = inline ('cos (x+y)');

dF1x = inline ('2*x+0*y');

dF1y = inline ('0*x+2*y');

disp (dF0x); disp (dF0y);

disp (dF1x); disp (dF1y);

% Начальное приближение

x=0; y=-1; e=0.1;

% Значение функций

root = newton (x, y, e, F0, F1,dF0x, dF0y, dF1x, dF1y);

disp ('Решение системы');

disp (root);

После выполнения программы получим следующее:

>>

Inline function:

F0(x, y) = sin (x+y)-1.6*x

Inline function:

F1(x, y) = x²+y²−1

Inline function:

dF0x (x, y) = cos (x+y)-1.6

Inline function:

dF0y (x, y) = cos (x+y)

Inline function:

dF1x (x, y) = 2*x+0*y

Inline function:

dF1y (x, y) = 0*x+2*y

Количество итераций

Решение системы

— 0.6163

— 0.7875

Полученное решение совпадает с рассчитанным.

Список используемой литературы

1. Н. Бахвалов, Н. Жидков, Г. Кобельков. Численные методы. М., 2002, 632 с.

2. Н. Калиткин. Численные методы. М., 1972,

3. А. Самарский.

Введение

в численные методы. М.,, 270с.

4. М. Лапчик, М. Рагулина, Е. Хеннер. Численные методы. М., 2004, 384с.

5. В. Потемкин. Система MATLAB. Справочное пособие. М., 1997, 350с.

6. Е. Алексеев, О. Чеснокова. MATLAB 7. М., 2006, 464с.