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

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

На заметку

Производная от полинома также является полиномом. Так, если некоторый

степени на единицу меньшей, чем исходный. Если коэффициенты исходного полинома формируют список a = [a1,a2,...,an,an +1 ], то коэффициенты полинома-производной формируют список b = [b1,b2,...,bn−1,bn ], причем имеет место соотношение bk = (n −k +1)ak для всех индексов k = 1,2,...,n . Встроенная функция polyder() по списку коэффициентов a исходного полинома вычисляет список b коэффициентов полинома-производной.

Командой dz=-polyval(a,z)/polyval(b,z) вычисляется приращение для переменной z. Это взятое со знаком минус отношение функции в текущей точке z к значению производной в этой же точке. Для вычисления значения полинома в точке мы использовали встроенную функцию polyval(). Первым аргументом функции передается список коэффициентов полинома. Второй аргумент функции – точка, в которой вычисляется значение полинома.

На заметку

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

После этого запускается оператор цикла. Проверяется условие abs(dz)>eps, которое состоит в том, что добавка dz (по абсолютной величине) к текущему значению переменной z превышает погрешность eps.

На заметку

Строго говоря, это не означает, что корень вычислен с нужной точностью. Другим словами, переменная eps в общем случае не определяет точность вычисления корня уравнения. Более того, такой алгоритм не обеспечивает гарантированного завершения итерационного процесса – если начальное приближение для корня выбрать неудачно, можем получить бесконечный цикл. Но такую экзотику здесь исследовать не будем.

В операторе цикла командой z=z+dz вычисляется новое приближение для корня, а затем командой dz=-polyval(a,z)/polyval(b,z) на основе этого нового приближения вычисляется добавка для следующего итерационного шага. После завершения оператора цикла в переменную x записываем полученный результат (команда x=z). На рис. 5.16 показано окно редактора m-файлов с кодом созданной функции.

221

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

Рис. 5.16. Код функции EqSolve3() в окне редактора m-файлов

РезультатывычислениякорнейполиномаспомощьюфункцииEqSolve3() приведены в документе на рис. 5.17.

Рис. 5.17. Результат вычисления корней полинома с помощью функции EqSolve3()

Здесь командой a=[1,-15,59,-45] создается список для исходного полинома – это в терминах переменной x выражение x3 −15x2 + 59x − 45 .

222

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

У этого полинома три корня: x = 1 , x = 5 и x = 9 . В последнем убедиться несложно – достаточно вычислить полином в соответствующих точках. Так, в результате выполнения команды polyval(a,[1,5,9]) получаем список из трех нулей.

Три корня полинома находим последовательно с помощью команд

EqSolve3(a,2,0.001), EqSolve3(a,6,0.001) и EqSolve3(a,12, 0.001). Результат, как видим, вполне приемлемый.

На заметку

Стоит напомнить, что в Matlab есть встроенная функция roots() для вычисления корней полинома. В этом отношении написание собственных функций для вычисления корней полиномов представляется малоперспективным.

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

function x=EqSolve4(f,x0,eps) z=x0;

h=eps/100;

n=1000; for i=1:n

df=(f(z+h)-f(z-h))/2/h; dz=f(z)/df;

z=z-dz; if(abs(dz)<eps) break; end

x=z; end

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

Переменной z присваивается в качестве значения начальное приближение для корня. Переменная h определяет шаг приращения аргумента для вычисления производной от функции уравнения. Значение этой переменной на два порядка меньше значения аргумента eps. Переменная n определяет максимально возможное количество итераций, выполняемых при поиске корня уравнения. Именно эта переменная указана в качестве верхней границы диапазона изменения индексной переменной.

В операторе цикла выполняется несколько простых команд. Командой df=(f(z+h)-f(z-h))/2/h вычисляется значение производной для функции в точке с аргументом z.

223

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

На заметку

Напомним, что по определению производной от функции f (z) называется предел f ′(z) = lim f (z +dz) − f (z) . В данном случае мы использовали

В принципе, чем меньше параметр h , тем точнее выражение для производной.

Командой dz=f(z)/df вычисляется величина декремента для корня уравнения. Командой z=z-dz вычисляется новое итерационное значение для корня уравнения. После этого выполняется условный оператор if(abs(dz)<eps) break. Суть этой команды в том, что если изменение значения корня на какой-то итерации становится меньше параметра eps, итерационный процесс прекращается. Таким образом, итерационный процесс продолжается до тех пор, пока изменение корня не станет достаточно малым или пока не будет выполнено n итераций.

Полученное в качестве итерационного процесса значение z возвращается результатом функции EqSolve4(). Окно редактора m-файлов с кодом функции представлено на рис. 5.18.

Рис. 5.18. Код функции EqSolve4() в окне редактора m-файлов

Примеры вызова функции EqSolve4() представлены в документе на рис. 5.19.

Командой f=@(x)(x^3-15*x^2+59*x-45) мы определяем функцию уравнения вида f (x) = x3 −15x2 + 59x − 45 . Это все тот же полином, у которого три корня: x = 1 , x = 5 и x = 9 . Находим все три корня этого уравнения командами EqSolve4(f,2,0.00001), EqSolve4(f,6,0.00001)

224

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

Рис. 5.19. Решение уравнения с помощью функции EqSolve4()

и EqSolve4(f,12,0.00001). С удовлетворением наблюдаем, что корни найдены.

Метод Ньютона, с некоторыми модификациями, может быть расширен для поиска решений систем алгебраических уравнений. Например, допустим, что нам нужно решить систему из n уравнений f1(x1,x2,...,xn ) = 0 ,

f2(x1,x2,...,xn ) = 0 , ..., fn(x1,x2,...,xn ) = 0 относительно переменных x1 , x2 , ..., xn . Для поиска решения системы необходимо знать начальные при-

ближения для каждой из переменных. На основе нулевых (начальных) приближений вычисляется первое приближение, затем второе, и так далее. Рассмотрим этот процесс подробнее. Допустим, что имеется k -е приближение для решения системы x1 = x1(k) , x2 = x2(k) , ..., xn = xn(k) . Для расчета следующего приближения необходимо вычислить матрицу производных:

Матрица вычисляется в точке k -го приближения. Другими словами, мы вы-

225