Следуйте за решением для уравнения, в то время как параметр изменяется

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

Параметрированное скалярное уравнение

Параметрированное уравнение, чтобы решить

sinh(x)-3x=a,

где a числовой параметр, который идет от 0 до 5. В a=0, одно решение этого уравнения x=0когда a не является слишком большим в абсолютном значении, уравнение имеет три решения. Чтобы визуализировать уравнение, создайте левую сторону уравнения как анонимная функция. Постройте функцию.

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)')

Figure contains an axes object. The axes object contains 2 objects of type line.

Когда a является слишком большим или слишком маленьким, существует только одно решение.

Основанный на проблеме Setup

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

x = optimvar('x');
expr = sinh(x) - 3*x;

Создайте и постройте решения

Запуск с начального решения x=0 в a=0, найдите решения для 100 значений a от 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')

Figure contains 2 axes objects. Axes object 1 contains an object of type line. Axes object 2 contains an object of type line.

Скачок в решении происходит рядом a=2.5. В той же точке, количестве вычислений функции, чтобы достигнуть решения увеличения от приблизительно 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

Figure contains an axes object. The axes object contains 31 objects of type line, text. These objects represent fun, Maximum a.

Как a увеличения, сначала решения перемещаются налево. Однако, когда a выше 2.5, около предыдущего решения больше нет решения. fzero требует, чтобы дополнительные вычисления функции искали решение и находит решение около x = 3. После этого значения решения перемещаются медленно направо как a увеличения далее. Решатель требует только приблизительно 10 вычислений функции для каждого последующего решения.

Выберите Different Solver

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')

Figure contains 2 axes objects. Axes object 1 contains 2 objects of type line. These objects represent Solution, Error of Solution. Axes object 2 contains an object of type line.

fsolve несколько более эффективно, чем fzero, требование приблизительно 7 или 8 вычислений функции на итерацию. Снова, когда решатель не находит решения около предыдущего значения, решатель требует, чтобы намного больше вычислений функции искало решение. На этот раз поиск неудачен. Последующие итерации требуют немногих вычислений функции по большей части, но не удаются найти решение. Error of Solution постройте показывает значение функции fun(x)-a.

Чтобы попытаться преодолеть локальный минимум не быть решением, ищите снова от различной стартовой точки когда 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')

Figure contains 2 axes objects. Axes object 1 contains an object of type line. Axes object 2 contains an object of type line.

На этот раз, fsolve восстанавливается с плохой начальной точки рядом a=2.5 и получает решение, похожее на то, полученное fzero. Количество вычислений функции для каждой итерации обычно равняется 8, увеличиваясь до приблизительно 30 в точке, где решение переходит.

Преобразуйте целевую функцию Используя fcn2optimexpr

Для некоторых целевых функций или версий программного обеспечения, необходимо преобразовать нелинейные функции в выражения оптимизации при помощи fcn2optimexpr. Смотрите Поддерживаемые Операции для Переменных и выражений Оптимизации и Преобразуйте Нелинейную Функцию в Выражение Оптимизации. В данном примере преобразуйте исходный функциональный fun используемый для графического вывода к выражению оптимизации expr:

expr = fcn2optimexpr(fun,x);

Остаток от примера является точно тем же самым после этого изменения в определении expr.

Смотрите также

| |

Похожие темы