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

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

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

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

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

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

x = optimvar('x');
expr = fcn2optimexpr(fun,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 функциональных оценок для каждого последующего решения.

Выберите 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')

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

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

| |

Похожие темы