Main Content

This example shows how to solve an equation repeatedly as a parameter changes by starting subsequent solutions from the previous solution point. Often, this process leads to efficient solutions. However, a solution can sometimes disappear, requiring a start from a new point or points.

The parameterized equation to solve is

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

where $$a$$ is a numeric parameter that goes from 0 to 5. At $$a=0$$, one solution to this equation is $$x=0$$. When $$a$$ is not too large in absolute value, the equation has three solutions. To visualize the equation, create the left side of the equation as an anonymous function. Plot the function.

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

When $$a$$ is too large or too small, there is only one solution.

To create an objective function for the problem-based approach, create an optimization expression `expr`

in an optimization variable `x`

.

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

Starting from the initial solution $$x=0$$ at $$a=0$$, find solutions for 100 values of $$a$$ from 0 through 5. Because `fun`

is a scalar nonlinear function, `solve`

calls `fzero`

as the solver.

Set up the problem object, options, and data structures for holding solution statistics.

prob = eqnproblem; options = optimset('Display','off'); sols = zeros(100,1); fevals = sols; as = linspace(0,5);

Solve the equation in a loop, starting from the first solution `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

Plot the solution as a function of the parameter `a`

and the number of function evaluations taken to reach the solution.

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 jump in the solution occurs near $$a=2.5$$. At the same point, the number of function evaluations to reach a solution increases from near 15 to near 40. To understand why, examine a more detailed plot of the function. Plot the function and every seventh solution point.

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

As $$a$$ increases, at first the solutions move to the left. However, when $$a$$ is above 2.5, there is no longer a solution near the previous solution. `fzero`

requires extra function evaluations to search for a solution, and finds a solution near `x = 3`

. After that, the solution values move slowly to the right as $$a$$ increases further. The solver requires only about 10 function evaluations for each subsequent solution.

The `fsolve`

solver can be more efficient than `fzero`

. However, `fsolve`

can become stuck in a local minimum and fail to solve the equation.

Set up the problem object, options, and data structures for holding solution statistics.

probfsolve = eqnproblem; sols = zeros(100,1); fevals = sols; infeas = sols; asfsolve = linspace(0,5);

Solve the equation in a loop, starting from the first solution `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

Plot the solution as a function of the parameter `a`

and the number of function evaluations taken to reach the solution.

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`

is somewhat more efficient than `fzero`

, requiring about 7 or 8 function evaluations per iteration. Again, when the solver finds no solution near the previous value, the solver requires many more function evaluations to search for a solution. This time, the search is unsuccessful. Subsequent iterations require few function evaluations for the most part, but fail to find a solution. The `Error of Solution`

plot shows the function value $$fun(x)-a$$.

To try to overcome a local minimum not being a solution, search again from a different start point when `fsolve`

returns with a negative exit flag. Set up the problem object, options, and data structures for holding solution statistics.

rng default % For reproducibility sols = zeros(100,1); fevals = sols; asfsolve = linspace(0,5);

Solve the equation in a loop, starting from the first solution `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

Plot the solution as a function of the parameter `a`

and the number of function evaluations taken to reach the solution.

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

This time, `fsolve`

recovers from the poor initial point near $$a=2.5$$ and obtains a solution similar to the one obtained by `fzero`

. The number of function evaluations for each iteration is typically 8, increasing to about 30 at the point where the solution jumps.

`fcn2optimexpr`

For some objective functions or software versions, you must convert nonlinear functions to optimization expressions by using `fcn2optimexpr`

. See Supported Operations on Optimization Variables and Expressions and Convert Nonlinear Function to Optimization Expression. For this example, convert the original function `fun`

used for plotting to the optimization expression `expr`

:

expr = fcn2optimexpr(fun,x);

The remainder of the example is exactly the same after this change to the definition of `expr`

.