MATLAB.Exponenta
–Û·Ë͇ Matlab&Toolboxes

Вычисления и приближение данных в MATLAB

Алгоритм,по которому работает функция fzero

Как уже говорилось в разделе Основные способы обращения к fzero, при обращении к fzero возможно указать:

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

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

Поиск в функции fzero интервала, содержащего корень

Если функция fzero вызвана с начальным приближением к корню

fzero(fun, x0)

а не с отрезком, на границах которого происходит смена знака исследуемой функции, то происходит поиск такого отрезка с целью отделения корня. Алгоритм отделения корня в функции fzero достаточно простой.

Сначала выбирается шаг равный (для не равного 0), или (для равного 0). Далее вычисляется значение исследуемой функции в точке и сравнивается со значением исследуемой функции в точке . Если и разных знаков, то корень отделен, а если знаки совпадают, то происходит расширение промежутка вправо от точка до и сравниваются знаки функции в точках a и b. Если они различны, то корень отделен, а если нет, то в цикле происходит последовательное увеличение шага в раз и расширение промежутка влево, а затем и вправо от точки с проверкой знаков функции на границах. Как только знаки станут различными, отделение корня прекращается и начинает работать алгоритм по его уточнению, описанный в разделе Алгоритм уточнения корня в fzero.

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

Такой алгоритм может и не отделить корень при неудачно заданном начальном приближении. Например, возьмем уравнение

(x-499)(x-498)=0

и зададим в качестве начального приближения значение =500. Проследим за каждым шагом алгоритма fzero. Для этого зададим в качестве третьего входного аргумента структуру opts с полем Display, значение которого установлено в 'iter'. Структуру opt сформируем при помощи функции optimset (см. раздел Входные аргументы fzero).

 
>> fun = inline('(x-499)*(x-498)');
>> opts = optimset('Display','iter');
>> x = fzero(fun, 500, opts)

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

Search for an interval around 500 containing a sign change:
 Func-count    a          f(a)             b          f(b)        Procedure
    1             500             2           500             2   initial interval
    3         485.858       159.574       514.142       244.426   search
    5             480           342           520           462   search
    7         471.716       717.147       528.284       886.853   search
    9             460          1482           540          1722   search
…

2031   -5.92547e+153  3.51112e+307  5.92547e+153  3.51112e+307   search
 2033   -8.37988e+153  7.02224e+307  8.37988e+153  7.02224e+307   search
 2035   -1.18509e+154  1.40445e+308  1.18509e+154  1.40445e+308   search
Exiting fzero: aborting search for an interval containing a s	ign change
    because NaN or Inf function value encountered during search.
(Function value at -1.67598e+154 is Inf.)
Check function or try again with a different starting value.
Задание же значения =498.5 приводит к удачному отделению одного из корней уже на втором шаге (равного 498, т.к. сначала происходит расширение промежутка влево) и его последующему уточнению.
>> x = fzero(fun, 498.5, opts)
 
Search for an interval around 498.5 containing a sign change:
 Func-count    a          f(a)             b          f(b)        Procedure
    1           498.5         -0.25         498.5         -0.25   initial interval
    2           484.4       198.552         498.5         -0.25   search
В разделе Основные способы обращения к fzero был рассмотрен пример уравнения

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

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

Алгоритм уточнения корня в функции fzero

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

  • метод половинного деления (bisection);
  • метод секущих (secant или linear interpolation);
  • метод обратной квадратичной интерполяции (inverse quadratic interpolation).

Приведем краткие сведения об этих методах.

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

  1. Задан отрезок [a,b] , на границах которого функция принимает значения разных знаков (предполагается, что функция непрерывна).
  2. Отрезок делится пополам точкой и из полученных половин [a,c] и [c,b] выбирается та, на границах которой знаки различны. Выбранная половина отрезка обозначается [a,b].
  3. Так продолжается до тех пор, пока длина полученного отрезка не станет меньше заданной точности . Поскольку корень непрерывной функции все время находится между точками a и b, то за приближенное с заданной точностью значение к корню можно взять либо точку a , либо точку b (см. рисунок).

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

Примечание.
Если для уточнения корня уравнения используются итерационные методы, которые по текущему приближению xk строят новое приближение xk+1 , то их можно сравнить по порядку сходимости. Под порядком сходимости итерационного метода понимают следующее. Если ek=xk-x* — ошибка на каждой итерации (здесь x* — точное решение уравнения), то порядком сходимости называют число p такое, что:

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

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

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

<
Выбор нового приближения xk+1 к корню по методу секущих

Итерационная формула метода секущих выглядит следующим образом

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

Теорема (о сходимости метода секущих).
Пусть

Тогда существует окрестность корня x* такая, что последовательность xk, построенная по методу секущих по любым двум различным точкам из этой окрестности, сходится к корню x* и существует предел

где

Рассмотрим теперь последний из трех численных методов, задействованных в алгоритме функции fzero.

Метод обратной квадратичной интерполяции обладает более высоким порядком сходимости, нежели метод секущих. Метод обратной квадратичной интерполяции на каждой итерации использует три предыдущих приближения к корню xk-2, xk-1 и xk для получения текущего приближения xk+1. Для этого строится квадратичная парабола

проходящая через точки

Абсцисса точки пересечения параболы X=G(y) c осью Ox выбирается в качестве следующего приближения xk+1 к корню уравнения. Поскольку построенная парабола имеет вид X=G(y), т.е. x является зависимой переменной, а y — независимой, то метод называется методом обратной квадратичной интерполяции.

Несложно получить формулы для коэффициентов a, b, и c этой параболы и далее определить точку пересечения с Ox (положив y=0), поскольку это обычный интерполяционный полином (см. раздел Интерполяционный полином в форме Лагранжа.), в котором независимая переменная есть y:

В методе обратной квадратичной интерполяции важно, чтобы значения , , были различны, иначе при вычислении нового приближения xk+1=G(0) произойдет деление на ноль.

Перейдем теперь к описанию алгоритма уточнения корня, реализованного в функции fzero. Он основан на алгоритме ZEROIN, который был разработан в Математическом центре Амстердама в 1960-х годах Ван Венгарденом и Деккером. Описание алгоритма приведено, например, в книге Дж. Форсайта, М. Мальколма и К. Моулера «Машинные методы математических вычислений», перевод которой на русский язык был опубликован в издательстве «Мир» в 1980г. В этой книге так же содержится фортрановская версия ZEROIN.

Приведем алгоритм fzero без некоторых вычислительных подробностей.

На каждом шаге итерационного алгоритма уточнения корня в fzero присутствуют три приближения a ,b и c к корню для поиска нового приближения (с самого начала значения a и c совпадают):

  • a — предыдущее приближение (равное b на предыдущей итерации);
  • b — последнее приближение;
  • c — предпоследнее приближение, как правило, лежащее по другую сторону от корня, чем точка b, т.е. такое, что знаки ƒ(b) иƒ(c) различны.

Если достигнута заданная точность (или значение ƒ(b) оказалось тождественно равным нулю), то происходит выход из цикла. Для проверки достижения заданной точности в функции fzero используется следующее условие

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

  • половинного деления;
  • секущих;
  • обратной квадратичной интерполяции.
При выборе метода принимается во внимание следующее. Если точки b и a близки, или , то следующее приближение к корню находится по методу половинного деления для точек c и b. Иначе проверяется, совпадают ли точки a и c . Если точки a и c совпадают, то следующее приближение к корню находится по методу секущих (по точкам b и a), а если точки a и c не совпадают, то следующее приближение вычисляется по трем точкам a, b и c по методу обратной квадратичной интерполяции. Если при этом получилась точка, лежащая не слишком близко к текущим, то она выбирается в качестве следующего приближения, а иначе берется точка, полученная половинным делением.

В следующем разделе мы рассмотрим принципиальное отличие функции fzero от функции fsolve (входящей в Optimization Toolbox), которая так же может быть применена для решения нелинейных уравнений.

Отличие fzero от fsolve для решения уравнений вида f(x)=0

Для решения систем нелинейных уравнений (и, в частности, одного уравнения) в Optimization Toolbox имеется функция fsolve. Существенным отличием fsolve от fzero является то, что

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

В качестве примера рассмотрим решение уравнения

которое, очевидно, не имеет корней, т.к для любого значения x. Создадим inline-функцию для вычисления левой части уравнения

>> fun = inline('cos(x)-exp(1e-3+x.^2)');

и применим fzero и fsolve с для определения корня.

Функция fzero не смогла отделить корень и выводит об этом сообщение в командное окно

>> x = fzero(fun, 0)
Exiting fzero: aborting search for an interval containing a sign change
    because NaN or Inf function value encountered during search.
(Function value at -28.9631 is -Inf.)
Check function or try again with a different starting value.
x =
   NaN

Функция fsolve находит нулевое значение

>> x = fsolve(fun, 0)
Optimization terminated: first-order optimality is less than options.TolFun.
x =
     0

Из графика видно, что именно в точке x=0 функция достигает своего наименьшего по абсолютной величине значения, не пересекая при этом ось абсцисс.


График функции функция .

Для получения такого графика в MATLAB с точной разметкой оси ординат приходится задать разметку так, как это сделано в приведенных ниже командах. Сначала координаты разметки, которые являются значением свойства YTick, записываются в вектор-строку YTicks. Для получения значения свойства используется функция get, ее входной аргумент gca возвращает указатель на текущие оси. Затем из вектор-строки YTicks транспонированием создается вектор-столбец и преобразуется в массив строк YLabels. Для преобразование вектора в массив строк применяется функция num2str, во втором ее входном аргументе задается формат представления чисел — c 10 цифрами после десятичной точки. Этот массив строк YLabels и задается в качестве разметки оси ординат. Если не предпринять таких действий, то все подписи к оси ординат будут одинаковы из-за недостаточного количества значащих цифр, выбираемого по умолчанию.

Команды для получения графика с точной разметкой

% создание графического окна белого цвета
figure('Color', 'w')
% построение графика функции на отрезке [-0.0001 0.0001]
fplot(fun, [-0.0001 0.0001])
% получение координат разметки
YTicks = get(gca, 'Ytick');
YTicks = YTicks';
% преобразование их в массив строк
YLabels = num2str(YTicks,'%12.10f');
% задание более точных координат разметки
set(gca, 'YTickLabel', YLabels)
% нанесение сетки
grid on
% вывод маркера на график, соответствующего значению функции в нуле
hold on
plot(0,fun(0), 'ro')

Поиск всех корней полиномов

Разумеется, вещественные корни полинома могут быть найдены при помощи функции fzero (см. раздел Основные способы обращения к fzero и следующие за ним). Однако, возникает вопрос об отделении корней (см. раздел Отделение корней ) и о том, как найти комплексные корни. В связи с этим, для нахождения всех (вещественных и комплексных) корней полиномов разработаны специальные методы, один из которых реализован в MATLAB в функции roots.

Работа функции roots основана на том, что корни полинома

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

В первую строку сопровождающей матрицы C размера n на n входят дроби, зависящие от коэффициентов полинома, а нижняя поддиагональ состоит из -1. Остальные элементы матрицы нулевые. Собственные числа сопровождающей матрицы C находятся в функции roots при помощи функции eig, так же входящей в состав вычислительных функций MATLAB.

Пример. Требуется найти корни полинома

Для этого создаем вектор p его коэффициентов (указывая все коэффициенты заданного полинома, включая нулевые):

>> p = [3 0 2 0 -1 4 1];

и вызываем функцию roots, которая возвращает вектор корней полинома

>> r = roots(p)

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

r =
      0.79356 +    0.62699i
      0.79356 -    0.62699i
     -0.20431 +     1.1899i
     -0.20431 -     1.1899i
     -0.94089              
     -0.23761  

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

>> res = polyval(p, r)

Видим, что полученные значения достаточно близки к нулю.

res =
  4.4409e-015 -9.5479e-015i
  4.4409e-015 +9.5479e-015i
 -5.9952e-015 -2.0151e-014i
 -5.9952e-015 +2.0151e-014i
    9.77e-015              
 -6.6613e-016  

Точные нули не получились из-за ошибок округления

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

Для подтверждения численной неустойчивости рассмотрим решение уравнения двенадцатой степени.

(x-1)(x-2)(x-3)(x-4)(x-5)(x-6)(x-7)(x-8)(x-9)(x-10)(x-10)(x-12)=0

Очевидно, что корни данного уравнения равны 1,2,K,12. Для использования функции roots нам нужно получить точные значения коэффициентов этого полинома, для этого обратимся к средствам Symbolic Math Toolbox. Зададим символьную переменную x и сам полином. Сами коэффициенты запишем в вектор p при помощи функции sym2poly и вызо вем функцию roots для нахождения всех корней полинома.

syms x
Q=(x-1)*(x-2)*(x-3)*(x-4)*(x-5)*(x-6)*(x-7)*(x-8)*(x-9)*(x-10)*(x-11)*(x-12)
p = sym2poly(Q)
format long
r = roots(p)

В результате получаем следующие корни

r =

  11.99999999566062
  11.00000002971861
   9.99999991236427
   9.00000014590805
   7.99999984893271
   7.00000010089053
   5.99999995644702
   5.00000001179920
   3.99999999813657
   3.00000000014623
   1.99999999999614
   1.00000000000000

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

Теперь изменим коэффициент при x11 в исходном полиноме на 10-5 и найдем корни полинома с такими коэффициентами.

q = p;
q(2) = q(2) + 1e-5;
r = roots(q)

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

r =
 11.67187782501375 + 0.20425035284474i
 11.67187782501375 - 0.20425035284474i
  9.48053399279692 + 0.74604616810422i
  9.48053399279692 - 0.74604616810422i
  7.36508954640540 + 0.21731250968556i
  7.36508954640540 - 0.21731250968556i
  5.96107641427432                    
  5.00408168821227                    
  3.99982673383837                    
  3.00000244088808                    
  1.99999999435455                    
  1.00000000000027

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

x = 0.99:0.01:12.01;
yp = polyval(p,x);
yq = polyval(q,x);
plot(x, yp, x, yq)
grid on
legend('exact polynom', 'polynom with error')


Поиск по сайту:

Система Orphus

Яндекс.Метрика