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

Параметризованное уравнение для решения:

$$\sinh (x)\mathrm{}-\mathrm{3x}$$= a,

где $$a$$ - числовой параметр, который находится в диапазоне от 0 до 5. При $$a\mathrm{=}$$0 одним решением этого уравнения $$\mathrm{является}\mathrm{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)')

Если $$\mathrm{значение\; слишком}$$ велико или слишком мало, существует только одно решение.

Установка на основе проблем

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

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

Создание и печать решений

Начиная с исходного решения $$x\mathrm{=}$$ 0 $$\mathrm{при}\mathrm{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\mathrm{=}\mathrm{}$$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

По мере $$$$увеличения, сначала решения перемещаются влево. Однако, когда $$а$$ выше 2,5, решение больше не существует вблизи предыдущего решения. fzero требует дополнительных оценок функций для поиска решения и находит решение вблизи x = 3. После этого значения раствора медленно перемещаются вправо $$\mathrm{по\; мере}$$дальнейшего увеличения. Решатель требует только около 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 график показывает значение функции $$\mathrm{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\mathrm{=}\mathrm{}$$2,5 и получает раствор, аналогичный полученному fzero. Количество оценок функций для каждой итерации обычно составляет 8, увеличиваясь примерно до 30 в точке скачка решения.

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

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

expr = fcn2optimexpr(fun,x);

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