Материал: Matlab. Практический подход. Самоучитель

Самоучитель Matlab

На заметку

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

С помощью команды hold on переходим в режим "удержания графика" – последующие графические построения будут выполняться в том же самом окне.

Затем запускается оператор цикла, в котором индексная переменная i пробегает значения от 1 до count. Эта индексная переменная нумерует итерации в вычислении корня. В первую очередь при выполнении итерации проверяется значение переменной state. При истинности значения выполняется группа команд, которыми формируются списки точек для отображения отрезков, выделяющих границы текущего интервала поиска решения, и хорды, соединяющей граничные точки на графике. Команды такие: L1y=[0,f(a)] (y - координаты отрезка, выделяющего левую границу интервала поиска решения), L1x=[a,a] (x - координаты отрезка, выделяющего левую границу интервала поиска решения), L2y=[0,f(b)] (y - координаты отрезка, выделяющего правую границу интервала поиска решения), L2x=[b,b] (x - координаты отрезка, выделяющего правую границу интервала поиска решения), Ly=[f(a),f(b)] (y - координаты границ хорды, соединяющей крайние точки на графике) и Lx=[a,b] (x - координаты границ хорды, соединяющей крайние точки на графике).

На заметку

Здесь следует учесть, что отрезки, которыми выделяется диапазон поиска решения, имеют такие координаты. Левый отрезок соединяет точки с координатами (a, 0) и (a, f (a)) , а правый отрезок соединяет точки (b, 0) и (b, f (b)). Хорда – это отрезок, который соединяет точки (a, f (a)) и (b, f (b)). Чтобы отобразить отрезок, достаточно указать его начальную и конечную точки. Технически это сводится к созданию двух списков с координатами начальной и конечной точек отрезка.

Командами plot(L1x,L1y,'LineWidth',2,'Color','blue'), plot (L2x,L2y,'LineWidth',2,'Color','blue') и plot(Lx,Ly,'Line Width',2,'Color','green') отображаются отрезки для выделения границ интервала поиска решения и хорда соответственно. Все линии отображаются с толщиной линии 2, вертикальные отрезки синим цветом, а хорда – зеленым. Напомним, что описанные выше построения выполняются в том случае, если значение аргумента state равно true. Следующая команда выполняется в любом случае: c=(a*f(b)-b*f(a))/(f(b)-f(a)). Этой командой вычисляется точка пересечения хорды с осью абсцисс. Затем необходимо в эту точку сместить одну из границ интервала поиска решения. Для этого в условном

216

Глава 5. Решение уравнений и оптимизация

операторе проверяем результат выражения f(a)*f(c)>0. Если результат истинен, то знак функции слева на границе интервала поиска решения и в "центре" (точке пересечения хорд с осью абсцисс) и сдвигается левая граница – для этого выполняем команду a=c. В противном случае сдвигается правая граница интервала поиска решения с помощью команды b=c. Блок оператора цикла по реализации итерационного процесса на этом завершается.

На заметку

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

В качестве результата функцией EqSolve2() возвращается значение, полученное как значение точки пересечения хорды с осью абсцисс на последней итерации. Поэтому выполняем команду x=c.

Следующая группа команд выполняется, если значение переменной state равно true. В частности, на графике на оси абсцисс выделяются точки начального интервала поиска решения и точка, найденная в качестве корня уравнения. Для отображения левой границы начального интервала поиска корня уравнения используем команду plot(intl(1),0,'ks','Marker FaceColor','k','MarkerSize',5. Здесь, поскольку в процессе поиска корня переменные a и b меняли свое значение, значение левой границы получаем инструкцией intl(1). Это x -координата отображаемой точки. Ее y -координата равна нулю. Инструкция 'ks' означает, что график (состоящий в данном случае из одной точки) отображается черным цветом и в качестве маркеров используются квадраты. Цвет заполнения маркера задается опцией 'MarkerFaceColor', значение 'k' означает черный цвет. Опцией 'MarkerSize' задается размер маркеров – в данном случае это 5.

Точка, найденная как корень уравнения, выделяется с помощью команды plot(x,0,'ko','MarkerFaceColor','k','MarkerSize',7). Используется маркер черного цвета круглой формы (инструкция 'ko'). Размер маркера установлен равным 7. Командой plot(intl(2),0,'ks', 'MarkerFaceColor','k','MarkerSize',5) правая граница интервала поиска решения выделяется так же, как и левая. Наконец, командой hold off отменяется режим удержания графика.

На рис. 5.13 показано командное окно с командами, которыми проверяется работоспособность созданного программного кода.

217

Самоучитель Matlab

Рис. 5.13. Командное окно с вызовом функции EqSolve2()

Рассмотрим команды подробнее. Сначала командой f=@(x)((x-1).* (x-5).*(x-9)) определяется функция, для которой ищется корень уравнения. В данном случае речь идет о функции f (x) = (x −1)(x − 5)(x − 9).

На заметку

Функция EqSolve2() описана так, что ее первый аргумент, являющийся ссылкой на функцию решаемого уравнения, должен применяться к аргументусписку. Такое предположение неявно подразумевалось при использовании команды plot(z,f(z),'LineWidth',3,'Color','red'), когда инструкцией f(z) вызывалась функция уравнения с аргументом-списком. Поэтому при определении функции уравнения мы использовали оператор поэлементного умножения *.

Командой EqSolve2(f,[0 4],100,false) определяется корень уравнения f (x) = 0 на интервале от 0 до 4 на основе 100 итераций без создания графика. В этом интервале находится корень x = 1 , который, собственно,

ивычисляется. Чтобы увидеть еще и графические построения, используем команду EqSolve2(f,[0 4],4,true). В данном случае корень вычисляем всего четырьмя итерациями – исключительно чтобы проследить последовательность построений. При большем числе итераций большинство вспомогательных линий будет на графике неразличимо. Помимо вычисления корня (с небольшой, откровенно говоря, точностью), создается еще

ирисунок, который представлен на рис. 5.14.

Метод, который рассмотрим далее, называется методом Ньютона или методом касательных. Он также применяется для решения уравнений вида f (x) = 0 и состоит в следующем.

Для поиска корня берется начальное значение для корня. В этой точке для графика функции f (x) строится касательная. Точка пересечения этой ка-

218

Глава 5. Решение уравнений и оптимизация

Рис. 5.14. Графика создана в результате вычисления корня методом хорд

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

Фактически речь идет о процессе, подразумевающем уточнение значения корня на основе значения, вычисленного на предыдущей итерации. Чтобы реализовать этот метод, необходимо знать не только функцию f (x), но и ее производную f ′(x). Кроме того, метод применим далеко не всегда.

На заметку

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

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

Если на n -м итерационном шаге получено приближение xn для корня уравнения, то на следующей итерации приближение для корня определя-

219

Самоучитель Matlab

касательных иллюстрируется на рис. 5.15.

Рис. 5.15. Решение уравнения методом касательных

Рассмотрим программный код функции, предназначенной для решения методом касательных уравнений вида f (x) = 0 , где функция f (x) является полиномом. Ниже приведен программный код этой функции:

function x=EqSolve3(a,x0,eps) z=x0;

b=polyder(a); dz=-polyval(a,z)/polyval(b,z); while(abs(dz)>eps)

z=z+dz; dz=-polyval(a,z)/polyval(b,z); end

x=z; end

Функция называется EqSolve3() и у нее три аргумента: массив a с коэффициентами исходного полинома, начальное приближение для корня x0 и параметр eps, имеющий отношение к точности вычисления корня. В качестве результата функцией возвращается корень уравнения (переменная x).

Локальной переменной z командой z=x0 в качестве значения присваивается начальное приближение для корня уравнения. Командой b=polyder(a) вычисляется список b, который представляет собой набор коэффициентов для производной от полинома, заданного списком a. Здесь мы использовали встроенную функцию polyder(), предназначенную для вычисления производной от полинома.

220