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

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

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

function x=NSolve(f,x0,N)

%Начальные приближения для поиска корней: z=x0';

%Шаг изменения аргумента для вычисления производной: h=0.00001;

%Количество независимых переменных:

n=length(x0);

%Единичная матрица: E=eye(n);

%Нулевой вектор-столбец: func=zeros(n,1);

%Нулевая матрица: A=zeros(n);

%Вложенные операторы цикла: for k=1:N

for i=1:n

F=f{i}; % указатель на функцию получаем из списка указателей func(i,1)=F(z'); % вычисление значения функции

for j=1:n

A(i,j)=(F(z'+h*E(:,j)')-F(z'))/h; % заполнение матрицы производных end

end

226

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

% Новое приближение для вычисления корня: z=z-A\func;

end

x=z'; % результат вычислений end

У функции NSolve() три аргумента:

•Список f с указателями на функции, которые задают систему алгебраических уравнений. Это особый список – список ячеек. Его мы будем обсуждать отдельно. Пока же для нас важно, что элементами списка f являются указатели на функции.

•Список x0 представляет собой набор с начальными значениями переменных, относительно которых решается система уравнений.

•Аргумент N задает количество итераций, которые выполняются при поиске решения уравнения.

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

Первой командой в функции мы объявляем локальную переменную и присваиваем ей значение: z=x0'. Переменная z получает в качестве значения транспонированный вектор (список) начальных значений x0. При присваивании значения вектор x0 транспонируется.

На заметку

При вычислениях мы используем векторы-столбцы – например, вектор значений независимых переменных z или вектор значений функций func (описывается далее). Вместе с тем, аргументами функции NSolve() передаются векторы-строки. Это делается для удобства ввода команд при вызове функции NSolve(). Поэтому периодически при вычислениях используется процедура транспонирования векторов.

Переменная h задает значение приращения аргумента при вычислении производной.

На заметку

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

Командой n=length(x0) в переменную n записывается размер списка начальных значений, то есть количество независимых переменных, относительно которых решается система уравнений. При этом неявно предполагается, что количество уравнений в системе такое же, как и количество независимых переменных.

227

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

Командой E=eye(n) создаем единичную матрицу, размер которой определяется количеством независимых переменных системы. Из этой единичной матрицы мы при вычислении частных производных будем извлекать векторы-столбцы (с единицей в соответствующей позиции). Также нам придется создать вектор-столбец со значениями функций системы уравнений. В качестве нулевого приближения создаем вектор-столбец, заполненный нулями. Достигается это командой func=zeros(n,1). Командой A=zeros(n) создается нулевая матрица – квадратная, размер задается количеством независимых переменных системы. После этого запускаются вложенные операторы цикла. Индексная переменная внешнего цикла пробегает значения от 1 до N и нумерует итерации вычисления корня. Первым внутренним циклом, индексная переменная которого пробегает значения от 1 до n (количество независимых переменных и характерный размер векторов и матриц, используемых в вычислениях), заполняется вектор-столбец со значениями функций системы и матрицы частных производных.

Во внутреннем операторе цикла командой F=f{i} переменной F в качестве значения присваивается указатель на функцию, индекс которой в списке указателей определяется значением индексной переменной i внутреннего цикла. Обратите внимание: вместо круглых скобок для индексирования элементов списка f используются фигурные скобки. Поступаем так, поскольку предполагаем (и не без оснований), что список f представляет собой список ячеек.

На заметку

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

С помощью команды func(i,1)=F(z') вычисляется значение соответствующей функции системы при заданных значениях независимых переменных. При этом аргумент, который представляет собой список значений независимых переменных, необходимо транспонировать. Полученный результат записывается в качестве соответствующего элемента списка func. Затем запускается еще один внутренний оператор цикла, в рамках которого заполняется матрица частных производных.

228

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

Основу тела цикла составляет команда A(i,j)=(F(z'+h*E(:,j)')-

F(z'))/h. Здесь для вычисления частной производной ∂F(x1,x2,...,xn )

∂xj

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

Чтобы решить первую задачу – вычислить значение функции с увеличенным на шаг дискретности одним аргументом, используем инструкцию F(z'+h*E(:,j)'). Здесь интерес представляет аргумент функции, который представляет собой сумму z'+h*E(:,j)'. Первое слагаемое тривиально – это транспонированный вектор текущих значений независимых переменных. Второе слагаемое – результат транспонирования вектора h*E(:,j). Инструкцией E(:,j) возвращается j-й столбец единичной матрицы E. У этого вектора-столбца все элементы нулевые, кроме j-го, который равен 1. После умножения на h получаем вектор-столбец со всеми нулевыми элементами, кроме j-го, который равен h. Поэтому в результате выполнения инструкции z'+h*E(:,j)' получаем вектор-столбец, у которого по сравнению с вектором z' элемент в j-й строке получает приращение h, чего мы и добивались. Прочие вычисления тривиальны.

Наконец, зная вектор значений функции и матрицы производных, можем вычислить следующее приближение для вектора значений независимых переменных. Этой цели служит команда z=z-A\func. После того, как все необходимые итерации выполнены, результат, после транспонирования, записываем в переменную x с помощью команды x=z'. На этом вычисление решения системы уравнений заканчивается. Окно редактора m-файлов с кодом созданной функции представлено на рис. 5.20.

Теперь воспользуемся созданной функцией для решения системы из двух уравнений: x2 + 2xy −x −14 = 0 и 3y + xy2 −27 = 0 . У этого уравнения есть два решения. Первое x = 2 и y = 3 , а также еще одно x ≈ 7.2353

иy ≈ −2.1502 . Эти решения и попытаемся найти. На рис. 5.21 представлено командное окно с примером создания списка указателей на функции

ивызовами функции NSolve() для решения соответствующей системы нелинейных алгебраических уравнений.

229

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

Рис. 5.20. Окно редактора m-файлов с кодом функции

Рис. 5.21. Командное окно с примерами использования функции NSolve() для решения системы нелинейных уравнений

230