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

Spline Toolbox

Решение уравнений, функции fzero и roots

Почему возникает задача о приближенном нахождении корней

Предположим, нам задано уравнение вида ƒ ( x ) = 0 , скажем

x2 + 3x − 1 = 0 .

Требуется найти корни уравнения, т.е. такие величины x* , что данное уравнение обращается в тождество, именно:

x2* + 3x* − 1 = 0 .

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

ax2 + bx + c = 0

по известным со школы формулам и находим ровно два корня (при обозначении корней уравнения иногда будем использовать индексы):

Корни могут и совпадать, если дискриминант D равен нулю, тогда будем говорить, что корень кратный (для квадратного уравнения кратности два). Обычно, когда говорят, что у уравнения сколько-то корней, то при этом учитывают их кратности.

В случае кубических уравнений ситуация не намного сложнее. Как было установлено в середине XVI века, уравнение приводится к уравнению , а его решение находится по формулам Кардано. Формулы для нахождения корней кубического уравнения, как и в случае квадратного уравнения, содержат радикалы, т.е. величины вида . Далее было показано, что уравнение 4-ой степени может быть сведено к решению квадратного и кубического уравнений. Однако, уравнения выше четвертой степени в общем случае не разрешимы в радикалах, что было доказано норвежским математиком Нильсом Хенриком Абелем в 1826г.

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

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

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

	figure
	axes
	fplot('1/x', [-20 -0.5])
	hold on
	fplot('1/x', [0.5 20])
	fplot('sin(x)', [-20 20], 'r')
	grid on


Графики функций sin x и .

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

В чем заключается задача о численном решении уравнений

В предыдущем разделе мы заметили, что в общем случае у уравнения вида f ( x ) = 0 может быть сколько угодно корней. Определим теперь, что значит приближенно найти какой-то корень. Обычно требуется найти корень с заданной точностью, обозначим ее ε . На практике это значение может быть каким угодно, т.е. например 10−3 , 10−11 и т.д. Предположим, что x* является точным корнем уравнения f ( x ) = 0 , т.е f ( x* ) = 0 . Разумеется, значение x* заранее неизвестно. Тогда задача о нахождении корня x* с точностью ε заключается в нахождении какой-либо точки , которая отстоит от точки x* не более чем на ε .

Сказанное демонстрируется на рисунке ниже. Для нас любая точка из отрезка будет корнем с точностью ε .


Задача о нахождении корня уравнения с заданной точностью.

Такая постановка задачи порождает сразу несколько вопросов. Первый из них — предположим, что в результате работы некоторого алгоритма мы нашли точку , как тогда определить ее расстояние до точного корня x* , значение которого заранее неизвестно. Достаточно маленькое значение f ( ) вообще говоря не значит, что точка близко к корню x* . Точно также, достаточно большое значение f ( ) вообще говоря не значит, что точка далеко от корня x* . Здесь в слова «достаточно маленькое значение» и «достаточно большое значение» не вкладывается точный смысл. Но, например, если , то это не вовсе не значит, что точка отстоит от x* не больше, чем на ε . Действительно, график функции может проходить близко к оси абсцисс вдали от корня. В этом убеждает пример, приведенный ниже на рисунке, где , но не является корнем с точностью ε .


Пример, в котором малость f ( ) ошибочно может быть принята за близость к x* .

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


Пример, в котором большое значение f ( ) ошибочно может быть принята за большое удаление от x* .

Второй вопрос — надо убедиться, что в ε-окрестности корня x* нет других корней. Иначе, если мы найдем , отстоящий от x* не более, чем на ε (и как-то проверим это), то будем ошибочно думать, что больше вблизи корней нет.

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

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

Отделение корней

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

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

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

Ясно, что данное уравнение не имеет корней, на отрезке [1, 2] его левая часть изменяет знак, но непрерывности у левой части уравнения на этом отрезке нет. Попытаемся решить данное уравнение при помощи функции fzero (про использование функции fzero написано в разделе Основные способы обращения к fzero и следующих за ним):

	x = fzero('1/(x-sqrt(2))', [1 2])
	x = 1.4142

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

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

Например, функция f (x) = x3 + x + 1 знак на границах отрезка [-1, 1], непрерывна на этом отрезке и монотонна, т.к. . Следовательно на отрезке [-1, 1] существует ровно один корень.

Рассмотрим, еще функцию . Она непрерывна на всей вещественной оси, но не монотонна. Действительно, ее производная меняет знак с минуса на плюс при переходе через точку x такую, что ex = 7 , т.е. при x = ln7 . Значит, при x < ln7 исследуемая функция убывает, а при x > ln7 она возрастает. Заметим, что , следовательно на отрезке есть ровно один корень. Аналогично, т.к. , то на отрезке есть так же ровно один корень. Итак, простое исследование убедило нас в наличии ровно двух корней у данной функции.

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

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

Предположим, что нам требуется отделить корни функции

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

	x = -5:0.1:5;
	y = x.^3+1/5*x.^2-212/25*x+204/25;

Для построения графика применим функцию plot, а для нанесения на него сетки — команду grid с аргументом on:

	plot(x, y)
	grid on

В нашем примере получается следующий график


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

Изучая поведение функции мы видим, что у нее есть три корня, первый около -3.5, а два других около 1.5. Однако нашей задачей было отделение корней, т.е. указание промежутков, на которых находится ровно один корень. Как видно из графика, для первого корня это промежуток [-4, -3]. Для того, чтобы отделить два других корня, удобно воспользоваться инструментом ZoomIn графического окна для увеличения масштаба просмотра вблизи точки x = 1.5 . После увеличения масштаба просмотра становится ясно, что промежутками, содержащими ровно один корень, являются отрезки [1, 1.4] и [1.8, 2.2]. Итак, все три корня отделены.

Следующая функция представляет более сложный пример для отделения корней

Ее график на отрезке [0, 3], построенный по равноотстоящим точкам с шагом 0.1, приведен ниже.


График функции , построенный с шагом 0.1

Из графика видно, что один корень принадлежит отрезку [1.5, 2.5]. А вблизи точки x = 1 ситуация не очень понятная, более того, увеличение масштаба ничего не дает. Выход состоит в построении графика вблизи точки x = 1 с более мелким шагом, причем в достаточно малой окрестности единицы. Например, график f ( x ) , построенный на отрезке [0.999, 1.002] с шагом 0.0001, дает хорошее представление о расположении двух оставшихся корней.


График функции , вблизи x = 1.

Становится понятным, что у данной функции есть два корня вблизи точки x = 1 . Один из них находится на отрезке [0.9995, 1.0005], а другой — на отрезке [1.0005, 1.0015]. Все три корня исследуемой функции отделены.

Заметим, что функция plot строит график табличной функции, заданной векторами x и y, которые мы получили по исходной непрерывной функции f ( x ) , соединяя точки отрезками. В зависимости от выбора точек для построения графика, он может оказаться слишком неточным или даже неверным. Приведем самый простой пример.

Предположим, требуется построить график функции

на отрезке [0.05, 1]. Сделаем это двумя способами. Первый из них состоит в вычислении функции с шагом 0.05 на данном отрезке и применении функции plot. Второй способ основан на использовании функции fplot. Для сравнения выведем результаты в одно графическое окно на разные оси, для чего применим функцию subplot

	x = 0.05:0.05:1;
	y = sin(1./x);
	figure
	subplot(2,1,1)
	plot(x,y)
	subplot(2,1,2)
	fplot('sin(1/x)', [0.05 1])


Графики, полученные при помощи plot (верхний) и fplot (нижний).

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

Рассмотрим еще один прием, полезный при отделении корней. Перед тем, как строить график исследуемой функции f ( x ) для отделения корней уравненияf ( x ) = 0 , иногда бывает полезно представить исходное уравнение f ( x ) = 0 в виде g ( x )=h ( x ) , где g ( x ) и h ( x ) — по возможности достаточно простые функции, поведение которых известно. Корни тогда находятся в абсциссах точек пересечения графиков этих функций.

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

Преобразуем его к эквивалентному виду

Ясно, что у данного уравнения есть только один корень — абсцисса точки пересечения графиков гиперболы и экспоненциальной функции ex . Отделить данный корень достаточно просто, например при x = 0.1 левая часть уравнения меньше нуля, а при x = 1 — больше нуля, следовательно на отрезке [0.1, 1] содержится ровно один корень.

Рассмотрим теперь уравнение

Совершив несложное преобразование (перенос единицы в правую часть и умножение обеих частей на xΠ ), приведем его к виду

На отрезке [-15, 15] строим графики функций . Находим точку пересечения, получается, что у исходного уравнения есть корень.

	figure
	fplot('sin(x)', [-15 15])
	hold on
	fplot('x-pi', [-15 15], 'r')
	grid on


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

Теперь строим график исходной функции на отрезке [-15, 15]. Он не пересекает ось абсцисс, следовательно корня нет.

	figure
	fplot('sin(x)/(x-pi)-1', [-15 15])
	grid on


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

Несложно понять, что в первом случае мы совершили неэквивалентное преобразование, которое привело к появлению корня x = Π у преобразованного уравнения. Поэтому при совершении перехода от уравнения вида f ( x )= 0 к уравнению вида g ( x )=h ( x ) требуется выполнять только эквивалентные преобразования для того, чтобы не изменить число корней.

Рассмотрим еще один пример — уравнение, зависящее от параметра:

Требуется провести отделение корней уравнения в зависимости от значений параметра a , изменяющихся от 0.2 до 1.3. Для этого удобно, во-первых, привести уравнение к виду . Во-вторых, можно построить графики левой и правой частей уравнения при изменении значения параметра a от 0.2 до 1.3 с некоторым шагом, скажем 0.1. Для этого удобно применить цикл for. В цикле меняется значение параметра a и строятся графики. В заголовке графика написано текущее значение . Команда pause означает, что надо нажать клавишу для продолжения работы программы (эти операторы удобно набрать в новом m-файле и выполнить, например выделив и нажав F9).

figure
x=-3:0.01:3;
for a = -0.2:0.1:1.3
    y1 = exp(-a*x);  
    y2 = a*x.^2; 
    plot(x, y1, x, y2)
    grid on
    title(['a=' num2str(a)])
    pause   
end

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


Различное число корней в зависимости от значения параметра .

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

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

Основные способы обращения к fzero

Для уточнения корня в MATLAB предназначена функция fzero. В самом простом варианте существует два способа обращения к ней:

  1. с указанием начального приближения к корню x0
    	 x = fzero(fun, x0)
  2. с указанием отрезка [a b], на котором отделен корень (см. предыдущий раздел Отделение корней)
    	x = fzero(fun, [a b])

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

Входной аргумент fun задает функцию. Он может быть строкой с выражением для функции, заключенной в апострофы, или inline-функцией, что удобно для достаточно простых функций. Кроме того, fun может быть:

  • указателем на функцию;
  • именем m-файла c функцией;
  • анонимной функцией;
  • вложенной функцией.

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

Рассмотрим простой пример. Требуется решить уравнение

cos x − x = 0

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

	format long e

(или в меню File выбираем пункт Preferences, в левой части окна Preferences переходим к разделу Command Window, а в правой части в списке Numeric format выбираем long e).

Теперь определим inline-функцию

	fun = inline('cos(x)-x');

и вызовем fzero для нахождения корня.

	x = fzero(fun, [0 pi/2])

Получаем

	x =
		7.390851332151607e-001

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

	fun(x-1e-16)
		ans =
	1.110223024625157e-016
	fun(x+1e-16)
		ans =
	-2.220446049250313e-016

Проверка показала, что знаки разные, следовательно корень найден с точностью 10−16 .

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

		x = fzero(fun, [1 2])
	??? Error using ==> fzero
	The function values at the interval endpoints must differ in sign.

Найдем теперь приближенное значение корня, указав функции fzero не отрезок, а начальное приближение к корню, например 1.

		x = fzero(fun, 1)
	x =
		7.390851332151607e-001

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

Рассмотрим другой пример, в котором требуется решить уравнение

В принципе ясно, что у данного уравнения корень только один, т.к. корень это абсцисса точки пересечения функций и h ( x ) = x2 − 10 . Построим графики функций g ( x ) и h ( x ) , например на отрезке [Π,3.5] :

	figure
	fplot('sqrt(x-pi)', [pi 3.5])
	hold on
	fplot('x^2-10', [pi 3.5], 'r')
	grid on


Графики функций

Из графиков видно, что абсцисса точки пересечения примерно равна 3.2. Выберем ее в качестве начального приближения для функции fzero:

	fun = inline('sqrt(x-pi)-x^2+10');
	x = fzero(fun, 3.2)

Exiting fzero: aborting search for an interval containing a sign change
    because complex function value encountered during search.
(Function value at 3.0609 is 0.63086+0.28406i.)
Check function or try again with a different starting value.
x =
   NaN

Мы обнаруживаем, что решение не найдено (NaN означает «не число» от Not a Number), а в командное окно вывелось предупреждение о прекращении поиска интервала, на границах которого функция изменяет знак, поскольку в процессе вычислений появилось комплексное значение. Это произошло потому, что функция fzero при поиске отрезка выбрала значение, не принадлежащее области определения нашей функции [Π,+∞] .

Правильным решением в данном примере является предварительное отделение корня. Из графиков g (x) и h (x) видно, что функция меняет знаки на [Π,1.25] . Укажем его в качестве второго аргумента fzero

	x = fzero(fun, [pi 3.25])

В этом случае получаем верный ответ

	x =
	    3.200386655543571e+000

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

Способы задания уравнения для fzero

В предыдущем разделе Основные способы обращения к fzero мы использовали только один способ для задания левой части уравнения f (x)=0 при помощи inline-функции. Такой способ удобен, если требуется найти корень достаточно простой функции. В этом разделе мы рассмотрим другие способы задания исследуемых функций.

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

 >> help elfun

Пример задания функции строкой.

>> x = fzero('x-exp(-x)', [0,1])
x =
  0.5671

Кроме того, исследуемая функция может быть задана при помощи анонимной функции (anonymous function). Разберем создание анонимных функций. В самом простом случае анонимная функция является функцией одного аргумента. Например, функция f(x)=x-e-x может быть задана при помощи анонимной функции следующим образом:

>> f = @(x) x-exp(-x)
f = 
    @(x) x-exp(-x)

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

Указатель на созданную анонимную функцию может быть теперь указан в качестве первого входного аргумента fzero:

>> x = fzero(f, [0,1])
x =
    0.5671

Анонимные функции позволяют искать корни функций, зависящих от параметров. Об этом написано в разделе Решение уравнения, зависящего от параметров.

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

Левая часть уравнения является функцией от x. Для вычисления этой функции при заданном x, т.е. для вычисления интеграла от 0 до x , воспользуемся функцией численного интегрирования quadl. Ее первым входным аргументом должна быть подынтегральная функция, которую можно задать строковым выражением. При записи выражения необходимо использовать поэлементные операции, в нашем примере — поэлементное умножение, задаваемое точкой со звездочкой. Вторым и третьим аргументами quadl являются, соответственно, нижний и верхний пределы интегрирования. Функцию quadl вместе со своими аргументами укажем в качестве выражения для анонимной функции:

 >> f = @(x)quadl('exp(-t).*sin(t)-exp(t).*cos(t)',0, x);

Для отделения корней построим график зависимости от x на отрезке [0, 3] при помощи fplot:

>> fplot(f,[0 3])
>> grid on

Из графика следует, что искомый корень заключен между 2 и 2.5. Уточним его при помощи функции fzero:

 >> x = fzero(f, [2 2.5])

Получаем результат:

x =
    2.2009

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

Правильное выражение для inline-функции таково:

 
>> f = inline('quadl(''exp(-t).*sin(t)-exp(t).*cos(t)'',0, x)')

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

Рассмотрим теперь примеры, в которых поиск корня связан с программированием функции.

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

В качестве примера рассмотрим следующую задачу, в которой требуется выяснить, при каком минимальное собственное число матрицы A=B3, где

станет равным 2. Запрограммируем функцию, которая по заданному x возвращает , где — минимальное собственное число матрицы A.Текст соответствующей функции eq3 приведен ниже (она должна быть сохранена в файле eq3.m, который находится в текущем каталоге или в каталоге, путь к которому известен MATLAB).

function f = eq3(x)
B = [x 1 1 1
       1 x 1 1
       1 1 x 1
       1 1 1 x];
A = B.^3;  
f = min(eig(A))-2;

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

>> fplot(@eq3, [0 2])
>> grid on

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

>> x = fzero(@eq3, [0 2])

Получаем результат:

x =
    1.4422

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

>> x = fzero('eq3', [0 2])

Результат будет тем же самым.

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

function x = eqsolve
x = fzero(@eq3, [0 2]);

function y = eq3(x);
y = 1-x;

После вызова основной функции

>> x = eqsolve
x =
     1	

получаем верное решение уравнения 1-x=0. Однако, если мы в fzero обратимся к функции не по указателю, а по имени, т.е. напишем

x = fzero(eq3, [0 2]);

вместо

x = fzero(@eq3, [0 2]);

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

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

Решение уравнения, зависящего от параметров

В этом разделе речь пойдет про решение уравнений вида

f(x;a1,a2,...,an)=0.

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

  • MATLAB: Mathematics: Function Functions: Parameterizing Functions Called by Function Functions: Providing Parameter Values to Anonymous Functions;
  • MATLAB: Mathematics: Function Functions: Parameterizing Functions Called by Function Functions: Providing Parameter Values Using Nested Functions;

Первый из них связан с анонимными функциями. В разделе Способы задания уравнения для fzero задействовали их для задания уравнений вида f(x)=0. Теперь рассмотрим, как анонимные функции могут быть применены для решения уравнений, в которые входят параметры.

Для примера возьмем простое уравнение

Анонимную функцию (зависящую только от x) можно определить так

f = @(x) exp(a1*x) - a2;

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

>> a1 = 1;
>> a2 = 2;
>> f = @(x) exp(a1*x) - a2;

приводит к корректному определению функции (в данном примере для a1=1 и a2=2). Заметим, что при a2 > 0 у данного уравнения есть только один корень. В принципе, теперь можно использовать функцию fzero для уточнения корня

>> x = fzero(f, 0)
x =
    0.6931

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

function x = solvepar(a1, a2)
f = @(x) exp(a1*x) - a2;
x = fzero(f, 0);

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

>> x = solvepar(1,2)
x =
    0.6931
>> x = solvepar(1.5,2)
x =
    0.4621

и т.д.

Еще один способ задания уравнения, который мы рассмотрим на примере предыдущего уравнения, состоит в использовании вложенных функций. Вложенная функция размещается внутри основной и имеет доступ к переменным основной функции (подробно про вложенные функции написано в разделе справочной системы MATLAB: Programming: Types of Functions: Nested Functions).

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

function x = solvenest(a1, a2)
    function f = fun(x)
        f = @(x) exp(a1*x) - a2;
    end
x = fzero(fun, 0);
end

При наличии вложенных функций обязательно завершать и вложенные и основную функции ключевым словом end. Во вложенной функции fun значения a1 и a2 известны, поскольку они являются входными аргументами основной функции solvenest. Вызов fzero происходит в основной функции, сама fzero обращается ко вложенной со значениями параметров, переданными основной функции solvenest. Для поиска корней достаточно вызвать функцию solvenest от значений параметров:

>> x = solvenest(1,2)
x =
    0.6931
>> x = solvenest(1.5,2)
x =
    0.4621

и т.д.

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

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

function f = myfun(x, a1, a2)
f = exp(a1*x) - a2;

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

>> x = fzero(@myfun, 0, [], 1, 2)
x =
    0.6931

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

Интерфейс fzero, получение информации о вычислениях, параметры вычислительного процесса

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

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

  • строковое выражением, заключенное в апострофы;
  • inline-функция;
  • указатель на функцию;
  • имя m-файла c функцией;
  • анонимная функцией;
  • вложенная функцией.

Выходные аргументы fzero

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

>> [x,f] = fzero('exp(-x)-x', [0 1])
x =
    0.5671
f =
  1.1102e-016

Для получения информации о завершении работы функции fzero ее следует вызывать с тремя выходными аргументами

[x, f, exitflag] = fzero(...)

При этом значения третьего аргумента exitflag соответствуют следующим ситуациям:

  • exitflag = 1 — итерационный процесс нахождения корня сошелся к корню;
  • exitflag = -1 — алгоритм был прерван заданной пользователем функцией, которая вызывается на каждой итерации в fzero

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

  • exitflag = -3 — при нахождении интервала, на границах которого исследуемая функция меняет знаки, значение функции оказалось Inf или NaN.
  • exitflag = -4 — при нахождении интервала, на границах которого исследуемая функция меняет знаки, значение функции оказалось комплексным.
  • exitflag = -5 — итерационный процесс сошелся, но возможно к точке, в которой функция имеет особенность.

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

В этом примере при отделении корня fzero вышла за область определения функции, т.к. подкоренное выражение стало меньше нуля и exitflag = -4. При этом в командное окно выводится еще и соответствующее сообщение.

>> fun = inline('sqrt(x-pi)-x^2+10');
>> [x, f, exitflag] = fzero(fun, 3.2)
x =
   NaN
f =
   NaN
exitflag =
    -4

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

>> [x, f, exitflag] = fzero('1/(x-sqrt(2))', [1 2])
x =
    1.4142
f =
  1.5012e+015
exitflag =
    -5

Функция fzero ищет только те корни, в которых график исследуемой функции пересекает ось абсцисс, а не касается ее (см. раздел Отличие fzero от fsolve. Если во входном аргументе fzero указан не отрезок, а начальное приближение, то fzero начинает отделение корня, расширяя отрезок. В приведенном ниже примере при отделении корня значение функции x2вышло за допустимые в типе double границы и процесс отделения корня был остановлен. Выходной аргумент exitflag получил соответствующее значение.

>> [x, f, exitflag] = fzero('x^2', 1)
x =
   NaN
f =
   NaN
exitflag =
    -3

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

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

[x, f, exitflag, output] = fzero(…)

В последний из них output записывается информация о ходе вычислений. Аргумент output является структурой со следующими полями:

  • output.algorithm — использованные алгоритмы;
  • output.funcCount — сколько раз была вычислена исследуемая функция;
  • output.intervaliterations — число итераций, сделанное для отделения корня;
  • output.iterations — число итераций, сделанных для уточнения корня;
  • output.message — информация о завершении алгоритма.

Например, при нахождении корня уравнения x-cos x=0 с начальным приближением 1 выводится следующая информация

>> [x, f, exitflag, output] = fzero('x-cos(x)', 1)
x =
    0.7391
f =
     0
exitflag =
     1
output = 
intervaliterations: 8
iterations: 5
funcCount: 21
algorithm: 'bisection, interpolation'
message: 'Zero found in the interval [0.68, 1.22627]'

Из этой информации, в частности, следует, что корень был отделен за 8 итераций, а для уточнения корня применялись методы bisection (половинного деления) и interpolation (обратной квадратичной интерполяции, или линейной интерполяции). Про реализованные в fzero методы и алгоритм нахождения корня написано в разделе Алгоритм, по которому работает fzero.

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

Входные аргументы fzero

До сих пор (см., например Способы задания уравнения для fzero) мы рассматривали способ вызова fzero с двумя входными аргументами:

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

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

  • Display — уровень отображаемой о ходе алгоритма информации. Значения: 'off' (не отображать), 'iter' (выводить информацию о каждой итерации), 'final' (выводить информацию о завершении), 'notify' (используется по умолчанию, предупреждать о получении расходящегося процесса).
  • FunValCheck — вывод сведений о значениях функции. Значения: 'on' (выводить предупреждение о получении комплексных значений или NaN), 'off' (по умолчанию сведения о значениях функции не выводятся).
  • OutputFcn — функция, вызываемая на каждом шаге алгоритма. Данную функцию следует запрограммировать отдельно (см. раздел Прерывание работы fzero, пример приложения с кнопкой Stop).
  • TolX — точность нахождения корня, которая оценивается по разности соседних приближений. Более подробно об этом написано в разделе Алгоритм, по которому работает fzero.

Для формирования структуры следует вызвать функцию optimset, указав в ее входных аргументах пары: НазваниеПоля, значение. Например, мы хотим найти корень уравнения x-cos x=0 с точностью 10 -4 . Для этого определяем inline-функцию, формируем структуру options (название структуры может быть и другим, например opt3) и вызываем fzero:

>> fun = inline('cos(x)-x');
>> options = optimset('TolX', 1e-4);
>> x = fzero(fun, 1, options)
x =
    0.7391

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

>> fun(x-1e-4)*fun(x+1e-4)
ans =
 -2.6884e-008

На самом деле способ проверки на достижение заданной точности, заложенный в алгоритме fzero, не всегда дает гарантированный результат. Например, если мы зададим точность 10-10 при решении уравнения x41=0, то мы не получим приближенное решение (точное здесь равно 0) с этой точностью:

>> options = optimset('TolX', 1e-10);
>> x = fzero('x^41', 1, options)
x =
 -1.1706e-008

Почему так происходит, разобрано в разделе Отличие fzero от fsolve.

В процессе нахождения корня иногда полезно проследить, что происходит на каждом шаге алгоритма. Для этого следует установить Display в 'iter'. Обратимся снова к рассмотренному выше и проследим за действиями fzero на каждом шаге:

>> fun = inline('sqrt(x-pi)-x^2+10');
>> options = optimset('Display', 'iter');
>> x = fzero(fun, 5, options)
 
Search for an interval around 5 containing a sign change:
 Func-count    a          f(a)             b          f(b)        Procedure
    1               5      -13.6368             5      -13.6368   initial interval
    3         4.85858      -12.2954       5.14142      -15.0201   search
    5             4.8      -11.7522           5.2      -15.6053   search
    7         4.71716      -10.9964       5.28284      -16.4451   search
    9             4.6      -9.95235           5.4      -17.6572   search
   11         4.43431      -8.52617       5.56569      -19.4199   search
   13             4.2      -6.61121           5.8      -22.0095   search
   15         3.86863      -4.11363       6.13137      -25.8646   search
   17             3.4      -1.05166           6.6      -31.7003   search
Exiting fzero: aborting search for an interval containing a sign change
    because complex function value encountered during search.
(Function value at 2.73726 is 2.5074+0.63587i.)
Check function or try again with a different starting value.

Мы видим, что fzero начала искать отрезок, на границах которого исследуемая функция изменяет знаки, отступая все больше влево и вправо от начального приближения 5 (в выводимой таблице границы обозначены a и b). Процесс отделения закончился неудачно, потому что значение 3.4 было уменьшено до 2.73726 при котором значение исследуемой функции стало комплексным (это видно из сообщения, выведенного перед завершением работы fzero).

Проанализируем также, почему в случае поиска корня уравнения x41=0 с точностью 10-10 на самом деле находится корень с меньшей точностью.

>> options = optimset('Display', 'iter', 'TolX', 1e-10);
>> x = fzero('x^41', 1, options)

Сначала fzero отделяет корень и успешно справляется с этой задачей за 24 вычисления исследуемой функции, в результате находится отрезок [-0.28, 1.9051], на границах которого функция принимает значения разных знаков:

Search for an interval around 1 containing a sign change:
 Func-count    a          f(a)             b          f(b)        Procedure
    1               1             1             1             1   initial interval
    3        0.971716      0.308396       1.02828       3.13791   search
  ….
   24           -0.28 -2.15516e-023        1.9051  2.99639e+011   search

Затем происходит уточнение корня

Search for a zero in the interval [-0.28, 1.9051]:
 Func-count    x          f(x)             Procedure
   24           -0.28 -2.15516e-023        initial
   25           -0.28 -2.15516e-023        interpolation
….
  102   -1.34016e-008  -1.4822e-323        bisection
  103    -1.1706e-008             0        interpolation

На 103-ей итерации значение функции оказывается точным нулем, т.к. вещественные числа типа double лежат в интервалах [-1.79769e+308, -2.22507e-308] и [2.22507e-308, 1.79769e+308]. Возведение же последнего полученного приближения -1.1706e-008 в 41-ую степень дает машинный ноль и вычисления прекращаются.

Примечание. Границы значений для вещественных типов double и single выводятся при помощи функций realmax, realmin:

>> realmax
ans =
  1.7977e+308
>> realmin
ans =
  2.2251e-308
>> realmax('single')
ans =
  3.4028e+038
>> realmin('single')
ans =
  1.1755e-038

Прерывание работы fzero, пример приложения с кнопкой Stop

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

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

Интерфейс управляющей функции следующий:

flag_stop = outfun(x, optimValues, state)

Функция fzero передает следующую информацию в управляющую функцию:

  • x — приближение к корню, полученное на текущей итерации;
  • optimValues — структура, содержащая информацию о текущей итерации;
  • state — текущий этап алгоритма.

Смысл полей структуры optimValues следующий

  • funccount — суммарное на данный момент количество вычислений исследуемой функции;
  • iteration — номер текущей итерации;
  • intervaliteration — число проделанных итераций для отделения корня;
  • fval — текущее значение исследуемой функции;
  • procedure — применяемая в данный момент процедура (отделение корня, интерполяция или половинное деление);
  • intervala —левая граница текущего интервала;
  • fvala — значение исследуемой функции на левой границе интервала;
  • intervalb —правая граница текущего интервала;
  • fvalb — значение исследуемой функции на правой границе интервала.

Значением входного аргумента state может быть одно из следующих:

  • ‘init’ — алгоритм fzero находится в стадии инициализации;
  • 'interrupt' — в fzero выполняется итерация;
  • 'iter' — в fzero происходит завершение итерации;
  • 'done' — итерации завершены, выполняется последний этап алгоритма fzero.

Приведем пример приложения с графическим интерфейсом пользователя solgui, в котором кнопка Stop служит для прерывания работы fzero. Окно приложения приведено на рисунке ниже. Оно содержит строку ввода для ввода исследуемой функции или имени вычисляющей ее m-функции, кнопку Start для начала вычислений и область вывода результата.

В основной функции приложения создается графическое окно, строка ввода, две кнопки и текстовая область. С событием Callback кнопки Start связывается функция btnStart_Callback, а событием Callback кнопки Stop — функция btnStop_Callback (программирование приложений с графическим интерфейсом пользователя рассматривается в разделе Приложения с GUI и дескрипторная графика

В функции btnStart_Callback происходит считывание содержимого строки ввода (т.е. функции), далее при помощи optimset формируется структура options, поле которой OutputFcn принимает значение указателя на управляющую функцию stop_by_user. Переменной flag_stop присваивается значение 0 и она сохраняется в данных приложения при помощи функции guidata. После этого вызывается функция fzero с третьим входным аргументом options. Таким образом, fzero знает, что stop_by_user является управляющей функцией. В зависимости от значения третьего выходного аргумента exitflag функции fzero (см. раздел Выходные аргументы fzero) либо выводится ответ, либо сообщение о том, что при нахождении корня возникли проблемы.

Управляющая функция stop_by_user реализована как вложенная функция (nested function) в функцию btnStart_Callback. В ней выполняется единственный оператор — из данных приложения извлекается значение флага остановки, которое и становится значением выходного аргумента управляющей функции stop_by_user. Единственное место в приложении, где значение флага изменяется — функция btnStop_Callback обработки события Callback кнопки Stop. В этой функции флагу присваивается единица и его значение сохраняется в данных приложения при помощи функции guidata.

function solgui
% основная функция приложения

% создание графического окна
figure('MenuBar', 'none',...
    'NumberTitle', 'off', ...
    'Name', 'solgui', ...
    'Position', [300 300 400 100])
% создание области ввода
uicontrol('Style', 'edit',...
    'position',[10 70 210 20],...
    'Tag', 'edtFun')
% создание кнопки Start
uicontrol('Style', 'pushbutton',...
    'position',[240 70 70 20],...
    'String', 'Start',...
    'Callback', @btnStart_Callback)
% создание кнопки Stop
uicontrol('Style', 'pushbutton',...
    'position',[320 70 70 20],...
    'String', 'Stop',...
    'Callback', @btnStop_Callback)
% создание области вывода результата
uicontrol('Style', 'text',...
    'position',[10 10 380 50],...
    'FontSize', 12,...
    'Tag', 'txtOuttext')
end

function btnStart_Callback(src, evt)
% функция обработки нажатия на кнопку Start
    function stop = stop_by_user(x, optimValues, state)
    % управляющая функция
        % чтение значения флага остановки
        stop =  guidata(src);
    end
    % получаем структуру указателей на объекты приложения
    handles = guihandles(src);
    % считываем в fun содержимое строки ввода
    fun = get(handles.edtFun, 'String');
    % формируем структуру options, заполняя поле с управляющей функцией
    options = optimset('OutputFcn', @stop_by_user);
    % пока присваиваем флагу остановки 0
    flag_stop = 0;
    % сохраняем значение флага остановки
    guidata(src, flag_stop)
    % вызываем fzero
    [x,f, exitflag] = fzero(fun, 1, options);
    % в зависимости от значения exitflag выводим результат или сообщение о том,
    % что корень не найден
    if exitflag == 1
        outtext = char(['x = ' num2str(x)], ['f = ' num2str(f)]);
    else
        outtext = 'There were some problems when finding the root';
    end
    set(handles.txtOuttext, 'String', outtext)
end

function btnStop_Callback(src, evt)
% функция обработки нажатия на кнопку Stop
    % привсаиваем 1 флагу остановки
    flag_stop = 1;
    % сохраняем его значение в данных приложения
    guidata(src, flag_stop)
end

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

function y = myfun(x)
y = sin(x);
pause(4)

Итак, после запуска приложения, вода имени myfun в строку ввода и нажатия на кнопку Start начинаются вычисления, которые останавливаются нажатием на кнопку Stop.


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

Система Orphus

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