В этом примере показано, как неоднократно решать уравнение, когда параметр изменяется путем запуска последующих решений с предыдущей точки решения. Часто, этот процесс приводит к эффективным решениям. Однако решение может иногда исчезать, требуя запуска от новой точки или точек.
Параметрированное уравнение, чтобы решить
,
где числовой параметр, который идет от 0 до 5. В , одно решение этого уравнения когда не является слишком большим в абсолютном значении, уравнение имеет три решения. Чтобы визуализировать уравнение, создайте левую сторону уравнения как анонимная функция. Постройте функцию.
fun = @(x)sinh(x) - 3*x; t = linspace(-3.5,3.5); plot(t,fun(t),t,zeros(size(t)),'k-') xlabel('x') ylabel('fun(x)')
Когда является слишком большим или слишком маленьким, существует только одно решение.
Чтобы создать целевую функцию для подхода, основанного на проблеме, создайте выражение оптимизации expr
в переменной x
оптимизации.
x = optimvar('x');
expr = sinh(x) - 3*x;
Запуск с начального решения в , найдите решения для 100 значений от 0 до 5. Поскольку fun
скалярная нелинейная функция, solve
вызовы fzero
как решатель.
Настройте объект задачи, опции и структуры данных для содержания статистики решения.
prob = eqnproblem; options = optimset('Display','off'); sols = zeros(100,1); fevals = sols; as = linspace(0,5);
Решите уравнение в цикле, начинающем с первого решения sols(1) = 0
.
for i = 2:length(as) x0.x = sols(i-1); % Start from previous solution prob.Equations = expr == as(i); [sol,~,~,output] = solve(prob,x0,'Options',options); sols(i) = sol.x; fevals(i) = output.funcCount; end
Постройте решение в зависимости от параметра a
и количество вычислений функции, взятых, чтобы достигнуть решения.
subplot(2,1,1) plot(as,sols,'ko') xlabel 'a' ylabel('Solution(x)') subplot(2,1,2) plot(fevals,'k*') xlabel('Iteration Number') ylabel('Fevals')
Скачок в решении происходит рядом . В той же точке, количестве вычислений функции, чтобы достигнуть решения увеличения от приблизительно 15 до приблизительно 40. Чтобы понять почему, исследуйте более подробный график функции. Постройте функцию и каждую седьмую точку решения.
figure t = linspace(-3.5,3.5); plot(t,fun(t)); hold on plot([-3.5,min(sols)],[2.5,2.5],'k--') legend('fun','Maximum a','Location','north','autoupdate','off') for a0 = 7:7:100 plot(sols(a0),as(a0),'ro') if mod(a0,2) == 1 text(sols(a0) + 0.15,as(a0) + 0.15,num2str(a0/7)) else text(sols(a0) - 0.3,as(a0) + 0.05,num2str(a0/7)) end end plot(t,zeros(size(t)),'k-') hold off
Как увеличения, сначала решения перемещаются налево. Однако, когда выше 2.5, около предыдущего решения больше нет решения. fzero
требует, чтобы дополнительные вычисления функции искали решение и находит решение около x = 3
. После этого значения решения перемещаются медленно направо как увеличения далее. Решатель требует только приблизительно 10 вычислений функции для каждого последующего решения.
fsolve
решатель может быть более эффективным, чем fzero
. Однако fsolve
может стать всунутым локальный минимум и может не решить уравнение.
Настройте объект задачи, опции и структуры данных для содержания статистики решения.
probfsolve = eqnproblem; sols = zeros(100,1); fevals = sols; infeas = sols; asfsolve = linspace(0,5);
Решите уравнение в цикле, начинающем с первого решения sols(1) = 0
.
for i = 2:length(as) x0.x = sols(i-1); % Start from previous solution probfsolve.Equations = expr == asfsolve(i); [sol,fval,~,output] = solve(probfsolve,x0,'Options',options,'Solver','fsolve'); sols(i) = sol.x; fevals(i) = output.funcCount; infeas(i) = fval; end
Постройте решение в зависимости от параметра a
и количество вычислений функции, взятых, чтобы достигнуть решения.
subplot(2,1,1) plot(asfsolve,sols,'ko',asfsolve,infeas,'r-') xlabel 'a' legend('Solution','Error of Solution','Location','best') subplot(2,1,2) plot(fevals,'k*') xlabel('Iteration Number') ylabel('Fevals')
fsolve
несколько более эффективно, чем fzero
, требование приблизительно 7 или 8 вычислений функции на итерацию. Снова, когда решатель не находит решения около предыдущего значения, решатель требует, чтобы намного больше вычислений функции искало решение. На этот раз поиск неудачен. Последующие итерации требуют немногих вычислений функции по большей части, но не удаются найти решение. Error of Solution
постройте показывает значение функции .
Чтобы попытаться преодолеть локальный минимум не быть решением, ищите снова от различной стартовой точки когда fsolve
возвращается с отрицательным выходным флагом. Настройте объект задачи, опции и структуры данных для содержания статистики решения.
rng default % For reproducibility sols = zeros(100,1); fevals = sols; asfsolve = linspace(0,5);
Решите уравнение в цикле, начинающем с первого решения sols(1) = 0
.
for i = 2:length(as) x0.x = sols(i-1); % Start from previous solution probfsolve.Equations = expr == asfsolve(i); [sol,~,exitflag,output] = solve(probfsolve,x0,'Options',options,'Solver','fsolve'); while exitflag <= 0 % If fsolve fails to find a solution x0.x = 5*randn; % Try again from a new start point fevals(i) = fevals(i) + output.funcCount; [sol,~,exitflag,output] = solve(probfsolve,x0,'Options',options,'Solver','fsolve'); end sols(i) = sol.x; fevals(i) = fevals(i) + output.funcCount; end
Постройте решение в зависимости от параметра a
и количество вычислений функции, взятых, чтобы достигнуть решения.
subplot(2,1,1) plot(asfsolve,sols,'ko') xlabel 'a' ylabel('Solution(x)') subplot(2,1,2) plot(fevals,'k*') xlabel('Iteration Number') ylabel('Fevals')
На этот раз, fsolve
восстанавливается с плохой начальной точки рядом и получает решение, похожее на то, полученное fzero
. Количество вычислений функции для каждой итерации обычно равняется 8, увеличиваясь до приблизительно 30 в точке, где решение переходит.
fcn2optimexpr
Для некоторых целевых функций или версий программного обеспечения, необходимо преобразовать нелинейные функции в выражения оптимизации при помощи fcn2optimexpr
. Смотрите Поддерживаемые Операции для Переменных и выражений Оптимизации и Преобразуйте Нелинейную Функцию в Выражение Оптимизации. В данном примере преобразуйте исходный функциональный fun
используемый для графического вывода к выражению оптимизации expr
:
expr = fcn2optimexpr(fun,x);
Остаток от примера является точно тем же самым после этого изменения в определении expr
.