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

Параметризованное уравнение, которое нужно решить,

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

Когда $$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')

Переход в решении происходит около $$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

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

Выберите другой решатель

The 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 вычислений функции на итерацию. Снова, когда решатель не находит решения вблизи предыдущего значения, решатель требует намного больше вычислений функции для поиска решения. На этот раз поиск неудачен. Последующие итерации требуют нескольких вычислений функции по большей части, но не могут найти решение. The 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')

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

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

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

expr = fcn2optimexpr(fun,x);

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