Метод ньютона для решения нелинейных уравнений скайлаб

6 Нелинейные уравнения и системы в SCILAB

Леонид Баршов 5 лет назад Просмотров:

2 —>poly([1 0 ],’x’,’c’) 1 + x Листинг 6. Рассмотрим примеры символьных операций с полином ами: —>p1=poly([-1 ],’x’,’c’) p1 = x —>p=poly([3-7 ],’x’,’c’) p = 3-7x + x —>p1+p //Сложение — 5x + x —>p1-p //Вычитание 4 + 9x — x —>p1*p //Умножение x — 16x + 4x —>p1/p //Деление x —>p1^ //Возведение в степень 1-4x + 4x —>p^(-1) //Возведение в отрицательную степень x + x Листинг 6.3 Функция roots(p) предназначена для решения алгебраического уравнения. Здесь p это полином созданный функцией poly и представляющий собой левую часть уравнения P(x)=0. Решим несколько алгебраических уравнений. ЗАДАЧА 6.1. Найти корни полинома x 4-8x 3 +8x -1=0.

3 Для решения этой задачи необходимо задать полином р. Сделаем это при помощи функции poly, предварительно определив вектор коэффициентов V. Обратите внимание, что в уравнении отсутствует переменная x в первой степени, это означает, что соответствующий коэффициент равен нулю: —>V=[ ]; —>p=poly(v,’x’,’c’) p = x — 8x + x Листинг 6.4 Теперь найдем корни полинома: —>X=roots(p) X =! !! !! !! ! Листинг 6.5 Графическое решение задачи 1, показанное на рис. 6.1 позволяет убедится, что корни найдены верно. 3 Рисунок 6.1. Графическое решение задачи 6.1 ЗАДАЧА 6. Найти корни полинома x x +0.6x -1=0. Решение этой задачи аналогично решению предыдущей, разница заключается в способе вызова необходимых для этого функций: —>roots(poly([ ],’x’,’c’)) 1 Графическим решением уравнения f(x)=0 является точка пересечения линии f(x) с осью абсцисс.

4 4! !! i!! i! Листинг 6.6 Не трудно заметить, что полином имеет один действительный ( рис. 6.) и два комплексных корня. Рисунок 6.. Графическое решение задачи 6.. ЗАДАЧА 6.3 Найти решение уравнения y(x)=0, если y(x)=x 4-18x +6. Решение этой задачи представлено в листинге 6.7 и отличается от предыдущих лишь способом определения полинома. Графическое решение представлено на рис >x=poly(0,’x’); —>y=x^4-18*x^+.6; —>roots(y)! !! !! !! ! Листинг 6.7

5 5 Рисунок 6.3. Графическое решение задачи Трансценден тные уравнения Уравнение f(x)=0, в котором неизвестное входит в аргумент трансцендентных функций, называется трансцендентным уравнением. К трансцендентным уравнениям принадлежат показательные, логарифмические, тригонометрические. В общем случае аналитическое решение уравнения f(x)=0 можно найти только для узкого класса функций. Чаще всего приходится решать это уравнение численными методами. Численное решение нелинейного уравнения проводят в два этапа. В начале отделяют корни уравнения, то есть находят достаточно тесные промежутки, в которых содержится только один корень. Эти промежутки называют интервалами изоляции корня, определить их можно, изобразив график функции f(x) или любым другим методом. На втором этапе проводят уточнение отделенных корней, или, иначе говоря, находят корни с заданной точностью. Для решения трансцендентных уравнений в Scilab применяют функцию fsolve(x0,f) где x0 начальное приближение, f функция, описывающая левую часть уравнения y(x)=0. Рассмотрим применение этой функции на примерах. ЗАДАЧА Найти решение уравнения ( 1) x = 0 x. Определим интервал изоляции корня заданного уравнения. Воспользуемся графическим методом отделения корней. Если выражение, стоящее в правой части уравнения представить в виде разности двух функций f(x) g(x)=0, то абсцисса точки пересечения линий f(x) и g(x) корень данного уравнения. В нашем случае 3 3 f ( x) = ( x 1), g ( x) = x. На рис. 6.4 видно, что корень данного уравнения лежит в интервале [0;1]. Выберем ноль в качестве начального приближения, зададим функцию, описывающую уравнение и решим его: —>deff(‘[y]=f1(x)’,’y1=((x-1)^)^(1/3),y=(x^)^(1/3),y=y1-y’) —>fsolve(0,f1) 0.5 Листинг 6.8 Методы определения интервала изоляции корня основаны на следующем свойстве: если непрерывная функции f(x) на интервале [ a, b] поменяла знак, то есть f(a) f(b) 6 6 Рисунок 6.4. Графическое решение задачи 6.4 ЗАДАЧА 6.5 Найти корни уравнения f(x)=e x /5-(x-1). На рис. 6.5 видно, что график функции f(x) трижды пересекает ось абсцисс, то есть уравнение имеет три корня. Рисунок 6.5. Графическое решение задачи 6.5 Последовательно вызывая функцию fsolve с различными начальными приближениями, получим все решения заданного уравнения: —>deff(‘[y]=f(x)’,’y=exp(x)/5-*(x-1)^’) —>x(1)=fsolve(0,f);x()=fsolve(,f);x(3)=fsolve(5,f); —>x x =! !! !! ! Листинг 6.9, Кроме того начальные приближения можно задать в виде вектора и тогда функцию можно вызвать один раз: —>fsolve([0;;5],f)! !! !

7 7 Листинг 6.10! ! ЗАДАЧА 6.6 Вычислить корни уравнения sin( x)-0.4x=0 в диапазоне [ 5π;5π]. Решение задачи представлено в листинге >deff(‘[y]=fff(x)’,’y=-0.4+sin(x)’) —>V=[-5*%pi:%pi:5*%pi]; X=fsolve(V,fff); —>X //Множество решений X =! ! Листинг 6.11 ЗАДАЧА 6.7 Найти решение уравнения y(x)=0, если y(x)=x 5 -x Не трудно заметить, что заданное уравнение полином пятой степени, который имеет один действительный корень ( рис. 6.6) в интервале от до 1. Рисунок 6.6. Графическое решение задачи 6.7 Решим эту задачу при помощи функции fsolve: —>deff(‘[f]=y(x)’,’f=x^5-x^3+1′) —>X=fsolve(-,y) X = Листинг 6.1 Теперь применим функцию roots: —>roots(poly([ ],’x’,’c’))! i!! i!! i!! i!! ! Листинг 6.13

8 8 Как видим, заданное уравнение, кроме действительного корня ( листинг 6.1), имеет и мнимые ( листинг 6.13). Поэтому для отыскания всех корней полинома лучше использовать функцию roots. 6.3 Сис т емы уравнений Если заданы m уравнений с n неизвестными и требуется найти последовательность из n чисел, которые одновременно удовлетворяют каждому из m уравнений, то говорят о системе уравнений. Для решения систем уравнений в Scilab так же применяют функцию fsolve(x0,f). ЗАДАЧА 6.8 Решить систему уравнений: < x +y =1; x 3 -y=0>. Графическое решение системы ( рис. 6.7) показывает, что она имеет две пары корней. Рисунок 6.7. Графическое решение системы уравнений Окружность и гипербола пересекаются в точках [0.8;0.6] и [-0.8;-0.6]. Эти значения приблизительны. Для того чтобы уточнить их, применим функцию fsolve, предварительно определив систему с помощью файл функции : function [y]=fun(x) y(1)=x(1)^+x()^-1; y()=x(1)^3-x(); endfunction —>exec(‘c:\fun.sce’); disp(‘exec done’); exec done —>fsolve([ ],fun) >fsolve([ ],fun) Листинг 6.14 ЗАДАЧА 6.9 В данной задаче исследуется система из трех нелинейных уравнений с тремя неизвестными: function [y]=fun(x) y(1)=x(1)^+x()^+x(3)^-1 y()=*x(1)^+x()^-4*x(3)

9 9 y(3)=3*x(1)^-4*x()+x(3)^ endfunction —>exec(‘d:\scilab 3\fun’);disp(‘exec done’); exec done —>fsolve([ ],fun)//решение системы! ! Листинг 6.15

Нелинейные системы и уравнения

В более общем случае мы имеем не одно уравнение (1), а систему нелинейных уравнений $$ \begin \tag <2>f_i(x_1, x_2, \ldots, x_n) = 0, \quad i = 1, 2, \ldots n. \end $$ Обозначим через \( \mathbf = (x_1, x_2, \ldots, x_n) \) вектор неизвестных и определим вектор-функцию \( \mathbf(\mathbf) = (f_1(\mathbf), f_2(\mathbf), \ldots, f_n(\mathbf)) \). Тогда система (2) записывается в виде $$ \begin \tag <3>\mathbf(\mathbf) = 0. \end $$ Частным случаем (3) является уравнение (1) (\( n = 1 \)). Второй пример (3) — система линейных алгебраических уравнений, когда \( \mathbf (\mathbf) = A \mathbf — \mathbf \).

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

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

При итерационном решении уравнений (1), (3) задается некоторое начальное приближение, достаточно близкое к искомому решению \( x^* \). В одношаговых итерационных методах новое приближение \( x_ \) определяется по предыдущему приближению \( x_k \). Говорят, что итерационный метод сходится с линейной скоростью, если \( x_ — x^* = O(x_k — x^*) \) и итерационный метод имеет квадратичную сходимость, если \( x_ — x^* = O(x_k — x^*)^2 \).

В итерационном методе Ньютона (методе касательных) для нового приближения имеем $$ \begin \tag <4>x_ = x_k + \frac, \quad k = 0, 1, \ldots, \end $$

Вычисления по (4) проводятся до тех пор, пока \( f(x_k) \) не станет близким к нулю. Более точно, до тех пор, пока \( |f_(x_k)| > \varepsilon \), где \( \varepsilon \) — малая величина.

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

Чтобы найти корень уравнения \( x^2 = 9 \) необходимо реализовать функции

Данная функция хорошо работает для приведенного примера. Однако, в общем случае могут возникать некоторые ошибки, которые нужно отлавливать. Например: пусть нужно решить уравнение \( \tanh(x) = 0 \), точное решение которого \( x = 0 \). Если \( |x_0| \leq 1.08 \), то метод сходится за шесть итераций.

Теперь зададим \( x_0 \) близким к \( 1.09 \). Возникнет переполнение

Возникнет деление на ноль, так как для \( x_7 = -126055892892.66042 \) значение \( \tanh(x_7) \) при машинном округлении равно \( 1.0 \) и поэтому \( f^\prime(x_7) = 1 — \tanh(x_7)^2 \) становится равной нулю в знаменателе.

Проблема заключается в том, что при таком начальном приближении метод Ньютона расходится.

Еще один недостаток функции naive_Newton заключается в том, что функция f(x) вызывается в два раза больше, чем необходимо.

Учитывая выше сказанное реализуем функцию с учетом следующего:

1. обрабатывать деление на ноль
2. задавать максимальное число итераций в случае расходимости метода
3. убрать лишний вызов функции f(x)

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

При реализации метода Ньютона нужно знать аналитическое выражение для производной \( f^\prime(x) \). Python содержит пакет SymPy, который можно использовать для создания функции dfdx . Для нашей задачи это можно реализовать следующим образом:

Решение нелинейных систем

Идея метода Ньютона для приближенного решения системы (2) заключается в следующем: имея некоторое приближение \( \pmb^ <(k)>\), мы находим следующее приближение \( \pmb^ <(k+1)>\), аппроксимируя \( \pmb(\pmb^<(k+1)>) \) линейным оператором и решая систему линейных алгебраических уравнений. Аппроксимируем нелинейную задачу \( \pmb(\pmb^<(k+1)>) = 0 \) линейной $$ \begin \tag <5>\pmb(\pmb^<(k)>) + \pmb(\pmb^<(k)>)(\pmb^ <(k+1)>— \pmb^<(k)>) = 0, \end $$ где \( \pmb(\pmb^<(k)>) \) — матрица Якоби (якобиан): $$ \pmb<\nabla F>(\pmb^<(k)>) = \begin \frac<\partial f_1(\pmb^<(k)>)> <\partial x_1>& \frac<\partial f_1(\pmb^<(k)>)> <\partial x_2>& \ldots & \frac<\partial f_1(\pmb^<(k)>)> <\partial x_n>\\ \frac<\partial f_2(\pmb^<(k)>)> <\partial x_1>& \frac<\partial f_2(\pmb^<(k)>)> <\partial x_2>& \ldots & \frac<\partial f_2(\pmb^<(k)>)> <\partial x_n>\\ \vdots & \vdots & \ldots & \vdots \\ \frac<\partial f_n(\pmb^<(k)>)> <\partial x_1>& \frac<\partial f_n(\pmb^<(k)>)> <\partial x_2>& \ldots & \frac<\partial f_n(\pmb^<(k)>)> <\partial x_n>\\ \end $$ Уравнение (5) является линейной системой с матрицей коэффициентов \( \pmb \) и вектором правой части \( -\pmb(\pmb^<(k)>) \). Систему можно переписать в виде $$ \pmb(\pmb^<(k)>)\pmb <\delta>= — \pmb(\pmb^<(k)>), $$ где \( \pmb <\delta>= \pmb^ <(k+1)>— \pmb^ <(k)>\).

Таким образом, \( k \)-я итерация метода Ньютона состоит из двух стадий:

1. Решается система линейных уравнений (СЛАУ) \( \pmb(\pmb^<(k)>)\pmb <\delta>= -\pmb(\pmb^<(k)>) \) относительно \( \pmb <\delta>\).

2. Находится значение вектора на следующей итерации \( \pmb^ <(k+1)>= \pmb^ <(k)>+ \pmb <\delta>\).

Для решения СЛАУ можно использовать приближенные методы. Можно также использовать метод Гаусса. Пакет numpy содержит модуль linalg , основанный на известной библиотеке LAPACK, в которой реализованы методы линейной алгебры. Инструкция x = numpy.linalg.solve(A, b) решает систему \( Ax = b \) методом Гаусса, реализованным в библиотеке LAPACK.

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

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

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

Доброго времени суток! В этой статье мы поговорим о решении систем нелинейных алгебраических уравнений в Matlab. Вслед за решением нелинейных уравнений, переходим к их системам, рассмотрим несколько методов реализации в Matlab.

Общая информация

Итак, в прошлой статье мы рассмотрели нелинейные уравнения и теперь необходимо решить системы таких уравнений. Система представляет собой набор нелинейных уравнений (их может быть два или более), для которых иногда возможно найти решение, которое будет подходить ко всем уравнениям в системе.

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

Чтобы решить СНАУ, необходимо воспользоваться итеративными методами. Это методы, которые за определенное количество шагов получают решение с определенной точностью. Также очень важно при решении задать достаточно близкое начальное приближение, то есть такой набор переменных, которые близки к решению. Если решается система из 2 уравнений, то приближение находится с помощью построение графика двух функций.

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

Оператор Matlab для решения СНАУ

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

Решить систему нелинейных уравнений с точность 10 -2 :
cos(x-1) + y = 0.5
x-cos(y) = 3

Нам дана система из 2 нелинейных уравнений и сначала лучше всего построить график. Воспользуемся командой ezplot в Matlab, только не забудем преобразовать уравнения к стандартному виду, где правая часть равна 0:

Функция ezplot строит график, принимая символьную запись уравнения, а для задания цвета и толщины линии воспользуемся функцией set. Посмотрим на вывод:

Как видно из графика, есть одно пересечение функций — то есть одно единственное решение данной системы нелинейных уравнений. И, как было сказано, по графику найдем приближение. Возьмем его как (3.0, 1.0). Теперь найдем решение с его помощью:

Создадим функцию m-файлом fun.m и поместим туда следующий код:

Заметьте, что эта функция принимает вектор приближений и возвращает вектор значений функции. То есть, вместо x здесь x(1), а вместо y — x(2). Это необходимо, потому что fsolve работает с векторами, а не с отдельными переменными.

И наконец, допишем функцию fsolve к коду построения графика таким образом:

Таким образом у нас образуется два m-файла. Первый строит график и вызывает функцию fsolve, а второй необходим для расчета самих значений функций. Если вы что-то не поняли, то в конце статьи будут исходники.

И в конце, приведем результаты:

xr (это вектор решений) =
3.3559 1.2069

fr (это значения функций при таких xr, они должны быть близки к 0) =
1.0e-09 *
0.5420 0.6829

ex (параметр сходимости, если он равен 1, то все сошлось) =
1

И, как же без графика с ответом:

Метод простых итераций в Matlab для решения СНАУ

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

1. Если возможно, строим график.
2. Из каждого уравнения выражаем неизвестную переменную след. образом: из 1 уравнения выражаем x1, из второго — x2, и т.д.
3. Выбираем начальное приближение X0, например (3.0 1.0)
4. Рассчитываем значение x1, x2. xn, которые получили на шаге 2, подставив значения из приближения X0.
5. Проверяем условие сходимости, (X-X0) должно быть меньше точности
6. Если 5 пункт не выполнился, то повторяем 4 пункт.

И перейдем к практике, тут станет все понятнее.
Решить систему нелинейных уравнений методом простых итераций с точность 10 -2 :
cos(x-1) + y = 0.5
x-cos(y) = 3

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

x-cos(y) = 3
cos(x-1) + y = 0.5

Далее приведем код в Matlab:

В этой части мы выразили x1 и x2 (у нас это ‘x’ и ‘y’) и задали точность.

В этой части в цикле выполняются пункты 4-6. То есть итеративно меняются значения x и y, пока отличия от предыдущего значения не станет меньше заданной точности.

k = 10
x = 3.3587
y = 1.2088

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

Метод Ньютона в Matlab для решения СНАУ

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

1. Если возможно, строим график.
2. Выбираем начальное приближение X0, например (3.0 1.0)
3. Рассчитываем матрицу Якоби w, это матрица частных производных каждого уравнения, считаем ее определитель для X0.
4. Находим вектор приращений, который рассчитывается как dx = -w -1 * f(X0)
5. Находим вектор решения X = X0 + dx
6. Проверяем условие сходимости, (X-X0) должно быть меньше точности

Далее, решим тот же пример, что и в предыдущих пунктах. Его график мы уже строили и начальное приближение останется таким же.
Решить систему нелинейных уравнений методом Ньютона с точность 10 -2 :
cos(x-1) + y = 0.5
x-cos(y) = 3

Перейдем к коду:

Сначала зададим начальное приближение. Затем необходимо просчитать матрицу Якоби, то есть частные производные по всем переменным. Воспользуемся символьным дифференцированием в Matlab, а именно командой diff с использованием символьных переменных.

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

В этой части кода выполняем первую итерацию, чтобы получить вектор решения и сравнивать его с вектором начального приближения. Отметим, чтобы посчитать значение символьной функции в Matlab, необходимо воспользоваться функцией subs. Эта функция заменяет переменную на числовое значение. Затем функция double рассчитает это числовое значение.

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

В этой части кода выполняется цикл по расчету решения с заданной точностью. Еще раз отметим, что если в первой итерации до цикла были использованы функции diff, double и subs для вычисления производной в Matlab, то в самом цикле матрица якоби была явно задана этими частными производными. Это сделано, чтобы показать возможности среды Matlab.

За 3 итерации достигнут правильный результат. Также важно сказать, что иногда такие методы могут зацикливаться и не закончить расчеты. Чтобы такого не было, мы прописали проверку на количество итераций и запретили выполнение более 100 итераций.

Заключение

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