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

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

ются эти самые вычисления с выполнения команды x=(b+a)/2, которой переменной x в качестве значения присваивается среднее значение для границ интервала. Это соответствует точке в центре интервала поиска решения. Далее реализуется такой алгоритм.

Проверяется длина интервала, на котором ищется корень. Если длина интервала не превышает удвоенное значение погрешности, процесс вычисления корня можно заканчивать. Для этого в условном операторе проверяется условие abs(b-a)<2*eps и в случае его истинности командой return завершается выполнение кода функции.

На заметку

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

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

Влюбом случае длина интервала поиска решения уменьшается вдвое, и задача, фактически, сводится к исходной: нужно найти корень на интервале, только интервал теперь в два раза меньше. Поэтому смело используем команду x=EqSolve(f,[a b],eps). Именно с этой командой связан рекурсивный вызов функции.

Вкомандном окне вводим и выполняем следующие инструкции (выделены жирным шрифтом):

>>f=@(x)((x-1)*(x-5)*(x-9));

>>x=EqSolve(f,[0 4],0.0001)

x =

0.9999

>>x=EqSolve(f,[2 4],0.0001)

??? Error using ==> EqSolve at 5

Неверно указаны границы интервала поиска корня!

>>x=EqSolve(f,[2 5.1],0.0001)

x =

5.0000

211

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

>> x=EqSolve(f,[-100 15],0.0001) x =

9.0000

Командой f=@(x)((x-1)*(x-5)*(x-9)) мы создаем функцию-полином f (x) = (x −1)(x − 5)(x − 9) и решать будем, соответственно, уравнение f (x) = 0 . У уравнения три корня: x = 1 , x = 5 и x = 9 . Эти корни будем искать. Первый корень x = 1 с точностью 0.0001 находим командой x=EqSolve(f,[0 4],0.0001). Команда x=EqSolve(f,[2 4],0.0001)

дает пример того, как неправильно указывать интервал для поиска решения: на интервале от 2 до 4 корней нет.

Два других корня находим командами x=EqSolve(f,[2 5.1],0.0001) и x=EqSolve(f,[-100 15],0.0001). В последнем случае указан интервал, на котором находится все три корня. То, что в результате найден корень x = 9 , в известном смысле является случайностью. Гипотетически могли получить любой из трех корней – дело в том, какой корень будет "отсеян" при делении интервала поиска корня. Результаты вычислений показаны на рис. 5.11.

Рис. 5.11. Вычисление корня уравнения методом половинного деления

212

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

На заметку

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

Следующий метод решения алгебраических уравнений, который мы рассмотрим, называется методом хорд. Он сильно напоминает метод половинного деления и состоит в следующем. Для поиска решения уравнения f (x) = 0 необходимо указать начальный интервал поиска. Необходимое условие, как и в методе половинного деления, состоит в том, чтобы на границах интервала поиска решения функция f (x) принимала значения разных знаков. Если это так, то точки на графике функции f (x) на границах интервала поиска решения соединяются отрезком (хордой). Определяется точка пересечения этой хорды с осью абсцисс, и вычисляется значение функции в этой точке. В нее смещается та из границ интервала поиска решения, для которой имеет место совпадение знаков функции. Процесс продолжается до достижения нужной точности. Процесс графически проиллюстрирован на рис. 5.12.

Рис. 5.12. Решение уравнения методом хорд

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

function x=EqSolve2(f,intl,count,state) n=100;

a=intl(1);

b=intl(2);

if(state) l=b-a;

z=a-0.1*l:1.2*l/n:b+0.1*l;

213

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

plot(z,f(z),'LineWidth',3,'Color','red'); grid on;

title('Решение уравнения вида f(x)=0 методом хорд',...

'BackgroundColor','yellow','FontSize',12,'FontWeight','bold',...

'Color','red','EdgeColor','red'); hold on;

end

for i=1:count if(state) L1y=[0,f(a)]; L1x=[a,a]; L2y=[0,f(b)]; L2x=[b,b]; Ly=[f(a),f(b)]; Lx=[a,b];

plot(L1x,L1y,'LineWidth',2,'Color','blue');

plot(L2x,L2y,'LineWidth',2,'Color','blue');

plot(Lx,Ly,'LineWidth',2,'Color','green'); end

c=(a*f(b)-b*f(a))/(f(b)-f(a)); if(f(a)*f(c)>0) a=c;

else b=c; end

end x=c;

if(state)

plot(intl(1),0,'ks','MarkerFaceColor','k','MarkerSize',5);

plot(x,0,'ko','MarkerFaceColor','k','MarkerSize',7);

plot(intl(2),0,'ks','MarkerFaceColor','k','MarkerSize',5); hold off;

end end

Функция называется EqSolve2(), у нее четыре аргумента, и она возвращает в качестве результата корень уравнения (соответствующая переменная обозначена как x). Функция, которая определяет решаемое уравнение, передается в виде указателя (обозначен как f) первым аргументом функции EqSolve2(). Второй аргумент (обозначен intl) функции EqSolve2() представляет собой список из двух элементов – через него реализуется интервал поиска корня. Третьим аргументом (переменная count) функции EqSolve2() передается целое число итераций, выполняемых при поиске решения уравнения. Четвертый, логический аргумент state, позволяет переходить в режим, при котором помимо вычисления корня также отображается график функции уравнения (с некоторыми дополнительными геометрическими построениями).

214

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

На заметку

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

Командой n=100 вводится переменная, определяющая количество точек, на основании которых будут строиться кривые. Переменным a и b командами a=intl(1) и b=intl(2) присваиваются начальные значения границ интервала поиска решения. Затем запускается условный оператор, в котором в качестве условия проверяется аргумент state. Равенство этого аргумента значению true означает, что помимо вычисления корня будем выполнять еще и геометрические построения. В данном случае строится график функции, выполняются настройки этого графика, а также выделяются вертикальными отрезками границы интервала поиска решения и первая хорда, которая соединяет граничные точки на графике функции.

Командой l=b-a вычисляется длина интервала поиска решения. Затем командой z=a-0.1*l:(b-a+0.2*l)/n:b+0.1*l формируется список со значениями аргумента функции для уравнения. При этом диапазон изменения аргумента для отображения на графике на 20% больше длины интервала поиска решения. Поэтому начальная точка для диапазона изменения аргумента вычисляется как a-0.1*l, а конечная – как b+0.1*l. Шаг дискретности равен 1.2*l/n. График строится командой plot(z, f(z),'LineWidth',3,'Color','red'). Опции, указанные при вызове функции plot(), означают, что толщина линии равна 3 и она красного цвета. Координатная сетка отображается командой grid on. С помощью функции title() задается заголовок графика. Первым аргументом функции передается текст 'Решение уравнения вида f(x)=0 методом хорд', который и отображается в качестве заголовка. Значение 'yellow' опции 'BackgroundColor' означает, что заголовок будет отображаться в области желтого цвета. Опция 'FontSize' со значением 12 задает размер шрифта для отображения заголовка. О том, что шрифт будет жирным, свидетельствует значение 'bold' для опции 'FontWeight'. Опцией 'Color' задается цвет шрифта, а опция 'EdgeColor' определяет цвет рамки вокруг области, в которой отображается заголовок. Значения этих опций указаны с помощью ключевого слова 'red', что означает красный цвет. Таким образом, заголовок будет отображаться красным цветом, жирным шрифтом, на желтом фоне в красной рамке.

215