Этот пример демонстрирует, что 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 function, преобразуйте эту систему в функцию 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 Численное решение системы нелинейных уравненийИспользование fsolve для системы уравнений, преобразованной в функцию MATLAB. Используйте начальную точку x (i) = -1,9 для нечетных индексов, и x (i) = 2 для четных индексов. Задайте 'Display' на 'iter' чтобы увидеть прогресс решателя. Задайте 'Jacobian' на 'on' использовать якобиан, заданный в 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. Ключевое различие между 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 The MathWorks, Inc.