Математическое моделирование. Решение нелинейных уравнений

Решение нелинейных уравнений

 

Любое нелинейное уравнение можно свести к виду

                                                                             (1.1)

Существует несколько методов решения этого уравнения, и применение того или иного метода определяется видом функции . При рассмотрении различных методов воспользуемся практически очевидным утверждением:

Если  - непрерывная функция на отрезке ;  и  (значения функции на концах отрезка) противоположного знака, то существует , такой, что , то есть  - корень уравнения (1.1).

Следствие: Если  не меняет знак на отрезке , то - единственный.

Это утверждение вытекает из того, что функция  в этом случае имеет монотонную производную (рис. 2).

 

Рис.2

 

Метод 1/2 деления (метод бисекции, метод дихотомии)

Основная идея метода половинного деления основывается на приведенном выше утверждении. Отрезок, содержащий корень , делится пополам и выбирается та его половина, на концах которой  принимает значения разных знаков. Далее процесс повторяется до тех пор, пока отрезок не станет достаточно мал, чтобы заранее заданная точность вычисления корня была достигнута. Скорость убывания, очевидно, равна скорости сходимости геометрической прогрессии со знаменателем 1/2.

 

Метод простой итерации

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

                                                                                                          (1.5)

Рис.5

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

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

                                                                                        (1.6)

На рис.5 показана геометрическая интерпретация метода простой итерации. Корень уравнения (1.4) есть абсцисса точки пересечения графиков функций  и . Пусть - начальное приближение. Так как  -график прямой, проходящей под 45^о к осям координат, то по построению ясно,

что . Построенная таким образом итерационная последовательность сходится к корню  уравнения (1.5).

Сходимость

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

Пусть в некоторой окрестности корня  уравнения (1.5) функция  дифференцируема и удовлетворяет неравенству

                                          ,                                                           (1.7)

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

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

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

Отсюда легко выводится аналогичное условие сходимости для метода Ньютона. Сравнивая формулы (1.2) и (1.6), можно представить решение уравнения (1.1) методом Ньютона как решение уравнения (1.5) методом простой итерации с функцией . Дифференцируя эту функцию, получим, что . Таким образом, условие сходимости для метода Ньютона выглядит как .

 

Для реализации итерационных методов приближенного решения нелинейных уравнений остается лишь запрограммировать формулы (1.2), (1.3), (1.4) или (1.6) получения последовательных приближений. Рекомендуется использовать цикл repeat...until, причем выход из цикла осуществить при выполнении условия , где -заданная точность.

 

Метод Гаусса

Рассмотрим систему уравнений N-го порядка:

 , где .                                        (2.1)

В матричном виде эта система записывается в виде . Метод Гаусса изучался в курсе высшей математики, и напомним, что основная идея метода заключается в следующем.

1. Прямой ход метода Гаусса: исключая  из всех уравнений, кроме первого, получаем в первом столбце коэффициентов матрицы нули во всех строках, кроме первой. Затем исключаем  из всех уравнений, кроме первого и второго. Продолжая эти действия дальше, исключаем , ,... и приходим, наконец, к матрице треугольного вида.

2. Обратный ход метода Гаусса: последовательно исключаются   из (N-1)-го уравнения,  из (N-2)-го уравнения и т.д. Это приводит к диагональному виду матрицы коэффициентов и, соответственно к решению задачи.

Для составления программы, реализующей метод Гаусса, удобно пользоваться приведенными ниже формулами. Здесь вводятся вспомогательные матрица C размерности  и вектор Y, состоящий из N элементов.

Сначала для прямого хода в цикле для К от 1 до N необходимо вычислять

                                            (2.2)

Затем обратный ход метода позволяет получить решение в виде                                                                              (2.3)

 

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

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

Замечание 3. Количество арифметических действий в методе Гаусса растет с увеличением порядка системы .

Метод Зейделя

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

 , где ; .                 (2.7)

 

Отличия двух методов хорошо иллюстрирует следующий пример:

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

Возьмем в качестве начального приближения . Тогда

 ,   ,   ,...

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

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

 

Достаточные условия сходимости метода Зейделя формулируются следующей теоремой:

Если  для ,                                                               (2.8)

причем хотя бы для одного  неравенство строгое, то итерационная последовательность (2.7) сходится к решению системы (2.1) при любом начальном приближении.

Докажем эту теорему для системы двух уравнений с двумя неизвестными:                

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

                                                                          (*)

для предыдущего шага итерации эта формула имеет вид

                                                                              (**) Вычтем (**) из (*) и возьмем разность по абсолютной величине: . Аналогично можно получить, что
. Известно, что итерационная схема сходится, если  при . Это означает, что при выполнении условий теоремы  и схема сходится, что и требовалось доказать.

Линейная аппроксимация

Ставится задача найти коэффициенты аппроксимирующего многочлена, если предположить, что он имеет линейный вид: . Тогда минимизируемая функция - функция двух переменных , и для нахождения наилучших по принципу наименьших квадратов коэффициентов получаем систему двух линейных алгебраических уравнений с двумя неизвестными:

                          .                                          (3.6)

Эта система может быть записана в виде

                                          ,

 где , , , , .                  (3.7)

Таким образом, вычисляя по этим формулам значения коэффициентов системы и свободных членов, можно найти затем ее решение каким-либо известным методом. Остается лишь запрограммировать формулы (3.7) перед тем, как воспользоваться стандартной процедурой решения систем линейных уравнений, описанной раньше, например. В результате будут найдены значения коэффициентов  и , а график искомой функции будет представлять собой прямую и может выглядеть, например, как на рис.6,в.

Квадратичная аппроксимация

Аналогично решается задача нахождения коэффициентов аппроксимирующего многочлена 2-го порядка. В этом случае он имеет вид  и график представляет собой параболу. Минимизируемая функция - функция трех переменных . Повторяя те же выкладки, что и в предыдущем разделе, можно прийти к системе трех уравнений с тремя неизвестными (проделайте их самостоятельно):

                                          ,

 где , , , ,

       , , , .                                          (3.8)

Как видно из формул (3.8), полученная система симметрична и хорошо решается методом Гаусса, не требуя выбора главных элементов. Однако для многочленов более высокой степени () использование такого метода не рекомендуется, так как вычисленные на ЭВМ параметры могут оказаться полностью искаженными ошибками округления.

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

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

 

Численное интегрирование

 

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

 

Определение: квадратурной формулой называется приближенное равенство вида                 ,                                            (4.1)

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

 

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

 

Рис.7

 

Разобьем отрезок  на элементарные отрезки  точками . Интеграл разобьется при этом на сумму элементарных интегралов , где , что соответствует разбиению площади исходной криволинейной трапеции на сумму площадей элементарных криволинейных трапеций (рис.7). Введем обозначения: , где  - середина элементарного отрезка. Для простоты шаг  будем считать постоянным ().

Формулы прямоугольников

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

 

Рис.8

 

Если взять за высоту прямоугольника значение , то , тогда, суммируя все элементарные площади, получаем квадратурную формулу левых прямоугольников:

                                 .                                                           (4.2)

Во втором случае , и тогда приходим к квадратурной формуле правых прямоугольников:

                                 .                                                           (4.3)

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

                                 .                                                      (4.4)

Названия формул (4.2) и (4.3) следуют из геометрической иллюстрации (рис.8). В соответствии с этим приближенное равенство (4.4) иногда называют формулой центральных прямоугольников. Ниже, когда будут рассматриваться вопросы о точности квадратурных формул, выяснится, что она точнее на порядок по сравнению с (4.2) и (4.3). По этой причине для численного интегрирования гораздо чаще применяют (4.4) и называют ее просто формулой прямоугольников. Для всех трех полученных формул их применение означает приближенную замену площади исходной криволинейной трапеции площадью ступенчатой фигуры, состоящей из элементарных прямоугольников (рис.9 - для формулы центральных прямоугольников).

Рис.9

Формула трапеций

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

.                                              (4.5)

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

 

                

Рис. 10                                                     Рис. 11

 

Формула Симпсона (парабол)

Если площадь элементарной криволинейной трапеции заменить площадью фигуры, расположенной под параболой, проходящей через точки  ,  и , то получим приближенное равенство . Здесь  - интерполяционный многочлен 2-й степени с узлами . На рис.11 показана такая парабола, а график функции  на элементарном отрезке для большей наглядности изображен искривленным. Если построить интерполяционный многочлен по указанным трем точкам, он будет выражаться формулой

.

Интегрирование полученного выражения на элементарном отрезке приводит к равенству

.                                  

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

                        (4.6)

 

Замечание 1. Учитывая геометрическую интерпретацию формулы Симпсона, ее иногда называют формулой парабол.

Замечание 2. В случае, когда число элементарных отрезков разбиения четно (), в формуле Симпсона можно использовать узлы лишь с целыми индексами:

                                                             (4.7)

Эта формула легко вытекает из (4.6), если при выводе принять  за элементарный отрезок длины .

Оценка погрешности

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

Оценим погрешность выведенных квадратурных формул. Будем использовать обозначение   - максимальное значение производной -го порядка от подынтегральной функции на отрезке . Выведем сначала оценку погрешности для формулы (4.2). Представим ее остаточный член в виде:

                              .

Пусть функция  непрерывно дифференцируема. Тогда, используя формулу Тейлора , где , имеем , и тогда

                   .

Так как , то, заменяя , приходим к оценке погрешности:

                                          .                                                             (4.8)

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

Подобным же способом оцениваются погрешности других квадратурных формул. Справедливы следующие утверждения (без доказательства):

 

Пусть функция  дважды непрерывно дифференцируема на отрезке интегрирования. Тогда для составной квадратурной формулы прямоугольников (центральных) справедлива оценка погрешности:

                                          ,                                                           (4.9)

а для формулы трапеций:

                                          .                                                   (4.10)

Приведенные оценки показывают, что формулы (4.2) и (4.3) имеют лишь первый порядок точности относительно , поэтому они используются крайне редко. Предпочтительнее использование формул (4.4) и (4.5), которые имеют второй порядок точности. Однако вычисление интегралов по формуле Симпсона дает четвертый порядок точности относительно  (для функции  здесь требуется существование непрерывной производной четвертого порядка):

                                          .                                                  (4.11)

Замечание3. Из оценок (4.9)-(4.11) следует, что формулы прямоугольников и трапеций точны для многочленов первой степени, а формула Симпсона - для многочленов четвертой степени.

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

 

Метод Эйлера

Простейший дискретный аналог дифференциального уравнения (5.1) представляет собой уравнение

                     ,
приводящее к методу Эйлера. Вычисление очередного значения  осуществляется здесь по формуле

                     .                                                                    (5.3)

Геометрическая интерпретация метода Эйлера проиллюстрирована на рис.12. На каждом шаге вычислений строится касательная к интегральной кривой в точке , и за  принимается ордината точки ее пересечения с
вертикальной прямой . В результате неизвестная интегральная кривая заменяется ломаной линией, называемой ломаной Эйлера.

Рис. 12

Реализация численного алгоритма на компьютере предполагает построение алгоритма, позволяющего вычислить решение поставленной дискретной задачи Коши. В данном случае (и во всех ниже рассматриваемых случаях) вычисления производятся по явной формуле - значение сеточной функции зависит только от предыдущих значений, и программирование алгоритма не вызовет особых затруднений. Такие методы называют явными. В неявных методах правая часть уравнения может зависеть от , и на каждом шаге возникает необходимость решения некоторого уравнения относительно . Примером здесь может служить неявный метод Эйлера, когда вычисления производят по формуле .

Метод Эйлера с уточнением

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

                                                                  (5.4)

Метод приближенного решения задачи Коши по формулам (5.4) относится к методам класса прогноз-коррекция и иногда его называют методом Эйлера-Коши. Геометрическая иллюстрация этого метода представлена на рис.13. Если принять ломаную Эйлера в виде цепочки векторов, то значение сеточной функции по формулам (5.4) можно представить как полусумму двух последующих звеньев-векторов.

Рис.13

Как видно из рисунка, вычисленное значение гораздо находится ближе к графику искомой функции. Метод Эйлера-Коши имеет второй порядок точности, и доказательство этого факта можно найти в [1].

 

Методы Рунге-Кутты

Наиболее популярными среди классических явных одношаговых методов являются методы Рунге-Кутты. Рассмотренные выше методы можно рассматривать как частные случаи этого класса методов. Самым известным из методов Рунге-Кутты является метод четвертого порядка точности. Вычисления на каждом шаге здесь проводятся последовательно в четыре этапа. Сначала вычисляются вспомогательные коэффициенты, а затем они подставляются в основную формулу:

                                                                     (5.5)

Этот метод весьма прост и, как показывает практика, довольно эффективен в обычных расчетах.

Литература

1. Амосов А.А., Дубинский Ю.А., Копченова Н.В. Вычислительные методы для инженеров. - М., Высшая школа, 1994.

Приложение: Варианты контрольных работ.

Вариант № 1

 

1. Методом половинного деления и методом хорд определить один из корней уравнения     с точностью до e=0,0001

2. Методом Гаусса найти решение системы уравнений:

                  

3. Подобрать по принципу наименьших квадратов для заданных значений X и Y квадратичную функцию:

x	-2	-1	0	1	2
y	24,43	12,86	7,02	7,74	14,60

4. Методами трапеций и Симпсона вычислить интеграл и сделать проверку:

5. Используя метод Эйлера 1^ого порядка, составить таблицу решения дифференциального уравнения  с начальным условием   и шагом  на отрезке

 

Вариант № 2

 

1. Методом половинного деления и методом касательных определить один из корней уравнения  с точностью до e=0,0001

2. Методом простых итераций найти решение системы уравнений с точностью e=0,00001:

                  

3. Подобрать по принципу наименьших квадратов для заданных значений X и Y квадратичную функцию:

x	-3	-2	-1	0	3
y	5,24	8,42	9,51	8,21	-6,62

4. Методами модиф. прямоугольников и трапеций вычислить интеграл и сделать проверку:

5. Используя метод Эйлера-Коши, составить таблицу решения дифференциального уравнения  с начальным условием   и шагом  на отрезке

 

Вариант № 3

 

1. Методом половинного деления и методом простых итераций определить один из корней уравнения  с точностью до e=0,0001

2. Методом Зейделя найти решение системы уравнений с точностью e=0,00001:

                  

3. Подобрать по принципу наименьших квадратов для заданных значений X и Y квадратичную функцию: