Multiple solutions using fsolve

Hi, I've been trying to solve a set of nonlinear equations in many variables. I see that multiple solutions exist for these and fsolve does not output all of them. I'm attaching the code for a simple example to illustrate this.

How do i extract all possible solutions from this?

Many Thanks,

Chandra.

The set of equations in the following example have 2 sets of solutions and fsolve outputs just one of them and gives an exit flag 1. On solving these equations by hand, i found that the solution to the variable a3 is quadratic and has 2 solutions which leads to a set of multiple solutions for all other variables.

% clear

% clc

function [] = nonlinsolve()

%%Solution from fsolve

options = optimset('Display','Iter');    

[sol,fval,exitflag] = fsolve(@sol_nonlin,zeros(1,15),options)

%%The solution i am looking for

a0 = 0;a(1)=a0;a(2)=a0;a(3)=2+2*sqrt(2);a(4)=2.6955;a(5)=3.6131;a(6)=a0;
a(7)=a0;a(8)=a0;a(9)=a0;a(10)=a0;a(11)=a0;a(12)=a0;a(13)=a0;a(14)=a0;
a(15)=a0;

f = sol_nonlin(a)

function f = sol_nonlin(a)

f(1) =  -a(2)^2/4;

f(2) = a(1) + a(2) - a(2)*a(5);

f(3) = - a(5)^2 + 2*a(5) + a(3) - (3*a(2)*a(9))/2 + 1;

f(4) = a(7) + 3*a(9) - 3*a(5)*a(9);

f(5) = a(2) - (a(2)*a(3))/2;

f(6) = a(3) + 2*a(4) + 2*a(5) - a(3)*a(5) - a(2)*a(7);

f(7) = 2*a(6) + 2*a(7) + 3*a(9) - (3*a(14)*a(2))/2 - (3*a(3)*a(9))/2 ...
       - 2*a(5)*a(7);
f(8) = 2*a(10) + 3*a(14) - 3*a(14)*a(5) - 3*a(7)*a(9);

f(9) = - a(3)^2/4 + a(3) - (a(2)*a(6))/2 + 1;

f(10) = a(6) + 2*a(7) + 3*a(8) - a(10)*a(2) - a(3)*a(7) - a(5)*a(6);

f(11) = - a(7)^2 + 2*a(10) + 3*a(11) + 3*a(14) - (3*a(15)*a(2))/2 ...
        - (3*a(14)*a(3))/2 - 2*a(10)*a(5) - (3*a(6)*a(9))/2;

f(12) = 3*a(12) + 3*a(15) - 3*a(15)*a(5) - 3*a(14)*a(7) - 3*a(10)*a(9);

f(13) = a(6) - (a(11)*a(2))/2 - (a(3)*a(6))/2;

f(14) = 2*a(10) + a(11) - a(12)*a(2) - a(10)*a(3) - a(11)*a(5) - a(6)*a(7);

f(15) = 2*a(12) + 3*a(15) - (3*a(13)*a(2))/2 - (3*a(15)*a(3))/2 ...
        - 2*a(12)*a(5) - (3*a(14)*a(6))/2 - 2*a(10)*a(7) ...
        - (3*a(11)*a(9))/2;

f(16) = 3*a(13) - 3*a(10)*a(14) - 3*a(13)*a(5) - 3*a(15)*a(7) ...
        - 3*a(12)*a(9);

Answer 1

[bym]

using a starting vector of

xo = [0 0 5 3 4 0 0 0 0 0 0 0 0 0 0]
[sol,fval,exitflag] = fsolve(@sol_nonlin,xo,options)

converges to your desired solution

sol =

Columns 1 through 8

   -0.0000   -0.0000    4.8284    2.6955    3.6131    0.0000    0.0000    0.0000

Columns 9 through 15

   -0.0000   -0.0000   -0.0000    0.0000   -0.0000   -0.0000    0.0000

[Chandrasekar Sureshkumar]

Thx for your response. It would be useful if you could tell me how you chose the initial guess without knowing the actual solution? This is a very simple case of a problem i'm trying to solve and randomly guessing the (correct) initial values each time is an impossible task.

Answer 2

[Walter Roberson]

fsolve() is not suitable for finding multiple solutions, except through the mechanism of running multiple times with different x0 until the number of distinct solutions found has accumulated to the number desired.

fmincon() of the square of sol_nonlin() can be a better tool (using the usual trick that zero solving can be expressed as minimizing the square of the function); you would work that by partitioning the permitted range of one of the variables, forcing it to adjust all the other variables to try to fit that one variable in to the target range. Knowing where and how finely to partition might not be simple, though, not unless you are allowed to use some prior knowledge or there is a useful approximation you are permitted to use.

[Chandrasekar Sureshkumar]

Thx for your response. I tried fmincon and the code is pasted below. I added the constraint function whose value is always positive. The solution is still not what i expect. The exit flag is 5 and i tried playing around with TolCon, TolX and TolFun, but there is still no change in the solution.

function [] = nonlinsolve()

clc

options = optimset('Display','Iter');

x0 = zeros(1,15);

% x0(3) = 5;x0(4)=3;x0(5)=4;

% [sol,fval,exitflag] = fsolve(@sol_nonlin,zeros(1,15),options)

% [sol,fval,exitflag] = fminunc(@sol_nonlin,zeros(1,15),options)

%% fmincon

[g,h] = con(x0);

f1 = sol_nonlin(x0);

[x,fval,exitflag] = fmincon('sol_nonlin',x0,[],[],[],[],[],[],'con',options)

con(x)

% %% The solution i am looking for

a = zeros(1,15);

a(3:5) = [4.8284 2.6955 3.6131];

function f1 = sol_nonlin(a)

f(1) = -a(2)^2/4;

f(2) = a(1) + a(2) - a(2)*a(5);

f(3) = - a(5)^2 + 2*a(5) + a(3) - (3*a(2)*a(9))/2 + 1;

f(4) = a(7) + 3*a(9) - 3*a(5)*a(9);

f(5) = a(2) - (a(2)*a(3))/2;

f(6) = a(3) + 2*a(4) + 2*a(5) - a(3)*a(5) - a(2)*a(7);

f(7) = 2*a(6) + 2*a(7) + 3*a(9) - (3*a(14)*a(2))/2 - (3*a(3)*a(9))/2 ...

- 2*a(5)*a(7);

f(8) = 2*a(10) + 3*a(14) - 3*a(14)*a(5) - 3*a(7)*a(9);

f(9) = - a(3)^2/4 + a(3) - (a(2)*a(6))/2 + 1;

f(10) = a(6) + 2*a(7) + 3*a(8) - a(10)*a(2) - a(3)*a(7) - a(5)*a(6);

f(11) = - a(7)^2 + 2*a(10) + 3*a(11) + 3*a(14) - (3*a(15)*a(2))/2 ...

- (3*a(14)*a(3))/2 - 2*a(10)*a(5) - (3*a(6)*a(9))/2;

f(12) = 3*a(12) + 3*a(15) - 3*a(15)*a(5) - 3*a(14)*a(7) - 3*a(10)*a(9);

f(13) = a(6) - (a(11)*a(2))/2 - (a(3)*a(6))/2;

f(14) = 2*a(10) + a(11) - a(12)*a(2) - a(10)*a(3) - a(11)*a(5) - a(6)*a(7);

f(15) = 2*a(12) + 3*a(15) - (3*a(13)*a(2))/2 - (3*a(15)*a(3))/2 ...

- 2*a(12)*a(5) - (3*a(14)*a(6))/2 - 2*a(10)*a(7) ...

- (3*a(11)*a(9))/2;

f(16) = 3*a(13) - 3*a(10)*a(14) - 3*a(13)*a(5) - 3*a(15)*a(7) ...

- 3*a(12)*a(9);

f1 = sum(f.^2);

function [g,h] = con(x)

f1 = 0;

xe = [0;0];

xo = [1;1];

N = 3;

V = [ f1 x(2) x(5) x(9);

x(1) x(3) x(7) x(14);

x(4) x(6) x(10) x(15);

x(8) x(11) x(12) x(13)];

tmp1 = 0;

for i = 0:N

tmp2 = 0;

for j = 0:N

tmp2 = tmp2 + V(i+1,j+1)*(xo(1)-xe(1))^i*(xo(2)-xe(2))^j;

end

tmp1 = tmp1 + tmp2;

end

g = -tmp1; % Inequality Constraint

h = []; % Equality Constraint