Этот пример демонстрирует, что Symbolic Math Toolbox помогает минимизировать ошибки при решении нелинейной системы уравнений.
Этот пример использует следующую поддержку Symbolic Math Toolbox:
Вычисление аналитического градиента с помощью gradient
Вычисление аналитического якобиана с помощью jacobian
Преобразование символьных результатов в числовые функции с помощью matlabFunction
Аналитическая визуализация с fplot
Цель состоит в том, чтобы найти минимум Функции Розенброка или так называемую "банановую" функцию.
Используйте fplot
создать быструю визуализацию Функции Розенброка двух переменных
syms x y a=1; b=100; f(x,y)=(a-x)^2+b*(y-x^2)^2
f(x, y) =
fsurf(f,[-2 2 -1 3])
view([-45 50])
colormap jet
caxis([0,1000])
colorbar
“Нелинейные уравнения с Аналитическим якобиевским” примером в Optimization Toolbox решают ту же задачу при помощи fsolve
функция. Рабочий процесс, показанный в том примере, имеет два потенциальных источника ошибок. Вы должны были бы
Переведите уравнения для Функции Розенброка и якобиана от математической формы в тексте к цифровому коду.
Явным образом вычислите якобиан. Для сложных уравнений эта задача подвержена ошибкам.
Этот пример показывает, что использование символьной математики, чтобы описать проблемный оператор и последующие шаги минимизирует или даже устраняет такие ошибки.
Во-первых, преобразуйте Функцию Розенброка f
к системе нелинейных уравнений F
, где F
градиент f
.
clear
n = 64;
x = sym('x',[n,1]);
f = sum(100*(x(2:2:n)-x(1:2:n).^2).^2 + (1-x(1:2:n)).^2);
F = gradient(f,x);
Покажите первые 10 уравнений:
F(1:10)
ans =
Чтобы совпадать с промежуточными результатами, показанными в примере Optimization Toolbox, используйте ту же форму F
.
F(1:2:n) = simplify(F(1:2:n) + 2*x(1:2:n).*F(2:2:n)); F(1:2:n) = -F(1:2:n)/2; F(2:2:n) = F(2:2:n)/20;
Теперь системы уравнений идентичны:
F(1:10)
ans =
Используйте jacobian
вычислить якобиан F
. Эта функция вычисляет якобиан символически, таким образом избегая ошибок, сопоставленных числовыми приближениями производных.
JF = jacobian(F,x);
Покажите первые 10 строк и столбцов якобиевской матрицы.
JF(1:10,1:10)
ans =
Большинство элементов якобиевского матричного JF
нули. Тем не менее, когда вы преобразуете эту матрицу в функцию MATLAB, результатом является плотная числовая матрица. Часто, операции на разреженных матрицах более эффективны, чем те же операции на плотных матрицах.
Поэтому, чтобы создать эффективный код, преобразуйте только ненулевые элементы JF
к символьному вектору прежде, чем вызвать matlabFunction. is
и js
представляйте шаблон разреженности JF
.
[is,js] = find(JF); JF = JF(JF~=0);
Система уравнений F
, представление Функции Розенброка, символьная матрица, которая состоит из символьных выражений. Смочь решить его с fsolve
функция, преобразуйте эту систему в функцию MATLAB. На этом шаге удобно преобразовать оба F
и его якобиан, JF
, к одной основанной на файле функции MATLAB, FJFfun
. В принципе это позволяет F
и JF
к переменным повторного использования. В качестве альтернативы можно преобразовать их в отдельные функции MATLAB.
matlabFunction([F;JF],'var',x,'file','FJFfun');
FJFfun
принимает переменные как список, но fsolve
ожидает вектор переменных. Функциональный genericObj
преобразует вектор в список. Анонимная функция bananaObj
задан, чтобы зафиксировать значения аргументов для genericObj
. Отметьте использование шаблона разреженности в функциональном genericObj
. Далее обратите внимание что этот подход использования FJFfun
и genericObj
вместе не связывается к конкретному выражению, представленному F
и может быть использован как есть для любого F
.
bananaobj = @(y) genericObj(y,@FJFfun,is,js)
bananaobj = function_handle with value:
@(y)genericObj(y,@FJFfun,is,js)
Используйте fsolve
для системы уравнений, преобразованной в функцию MATLAB. Используйте начальную точку x (i) = –1.9 для нечетных индексов и x (i) = 2 для ровных индексов. Установите 'Display'
to'iter'
видеть прогресс решателя. Установите 'Jacobian'
к 'on'
использовать якобиан задало in bananaobj
.
x0(1:n,1) = -1.9; x0(2:2:n,1) = 2; options = optimoptions(@fsolve,'Display','iter','Jacobian','on'); [sol1,Fsol,exitflag,output,JAC] = fsolve(bananaobj,x0,options);
Norm of First-order Trust-region Iteration Func-count f(x) step optimality radius 0 1 8563.84 615 1 1 2 3093.71 1 329 1 2 3 225.104 2.5 34.8 2.5 3 4 212.48 6.25 34.1 6.25 4 5 212.48 6.25 34.1 6.25 5 6 212.48 1.5625 34.1 1.56 6 7 116.793 0.390625 5.79 0.391 7 8 109.947 0.390625 0.753 0.391 8 9 99.4696 0.976563 1.2 0.977 9 10 83.6416 2.44141 7.13 2.44 10 11 77.7663 2.44141 9.94 2.44 11 12 77.7663 2.44141 9.94 2.44 12 13 43.013 0.610352 1.38 0.61 13 14 36.4334 0.610352 1.58 0.61 14 15 34.1448 1.52588 6.71 1.53 15 16 18.0108 1.52588 4.91 1.53 16 17 8.48336 1.52588 3.74 1.53 17 18 3.74566 1.52588 3.58 1.53 18 19 1.46166 1.52588 3.32 1.53 19 20 0.29997 1.24265 1.94 1.53 20 21 0 0.0547695 0 1.53 Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient.
Система нелинейных уравнений, F
, может быть альтернативно решен с vpasolve
функция, которая является числовым решателем, доступным в Symbolic Math Toolbox. Основное отличие between vpasolve
и fsolve
тот fsolve
использование арифметика с двойной точностью, чтобы решить уравнения, тогда как vpasolve
арифметика переменной точности использования и поэтому, позволяет вам управлять точностью расчетов.
sol2 = vpasolve(F);
Если вы присваиваете решение системы уравнений к одной переменной, sol2
, затем vpasolve
возвращает решения в форме массива структур. Можно получить доступ к каждому решению (каждое поле массива структур):
sol2.x1
ans =
Также можно преобразовать sol2
к символьному вектору, и затем получают доступ к области значений решений. Например, отобразите 5-е к 25-м решениям:
for k=(1:64) sol2_array(k) = sol2.(char(x(k))); end sol2_array(5:25)
ans =
function [F,J] = genericObj(x,FJFfun,is,js) % genericObj takes as input, vector x, Function and Jacobian evaluation % function handle FJFfun, and sparsity pattern is,js and returns as output, % the numeric values of the Function and Jacobian: F and Jfunction % FJs(1:length(x)) is the numeric vector for the Function % FJs(length(x)+1:end) is the numeric vector for the non-zero entries of the Jacobian xcell = num2cell(x); FJs = FJFfun(xcell{:}); F = FJs(1:length(x)); J = sparse(is,js,FJs(length(x)+1:end)); end
Copyright © 2016 MathWorks, Inc.