Why do I get different outputs for fsolve when I use different initial point?

Dear experts,

I was trying to solve the function cascade shown below using the fsolve on the second script. The problems uses a guess value for calculation of the functions on the first code below. The gues value will also be used on the end solution of the second code. But when I use different initial conditions, I get different answers. Can you please help me?

1. Below I put my two code scripts (one to call functions and the other to run the solver)

function [sol] = find_CP_MSNP(Cp_guess, phi_start, Jv,Z, z_pol,phi,X_mem)
Cp_guess_MS = [Cp_guess 1- sum(Cp_guess)]; 
[y0, yp0]  = decic(@MS_AA_imp, 0, [phi_start 0], [0 0 0 0 0], [1 1 0 0 0], [],[],Cp_guess_MS, Jv, Z);
[zsol_MS, ysol_MS] = ode15i(@(z,x,xp)MS_AA_imp(z,x,xp ,Cp_guess_MS, Jv, Z), [0 z_pol], y0, yp0);
zz=[0.00001 0.00001 0.00001 0.000001];
z=[0.00001 0.00001 0.00001];
C_surf = ysol_MS(end,:);

% define anonymous function 'f' for the partition function which accepts C_surf as an input parameter
f = @(zz)partitionFirst(zz,C_surf(1:3),Z,phi,X_mem);
% solve partitionFirst function to get C_entrance indirectly via anonymous function 'f'
C_entrance = fsolve(f,zz);

% define anonymous function 'nernstPlanckAnonymous' for the NernstPlanck_Vol function which accepts C_entrance as a parameter
% solve NernstPlanck_Vol function to get C_end indirectly via anonymous function 'nernstPlanckAnonymous'
nernstPlanckAnonymous = @(z)NernstPlanck_Vol(z,C_entrance(1:3));
C_end = fsolve(nernstPlanckAnonymous,z); %concentration on the end of membrane

%%second partition to calculate the permeate concentrate
secondpartitionAnonymous = @(z)partitionsecond(z, C_end(1:3)); %%%replaceing Co from partitionsecond to C_end for the next calculation of Cp_cal
Cp_calc = fsolve(secondpartitionAnonymous,z);

sol = Cp_calc(1:3) - Cp_guess(1:3);
end

2. The second code to solve this function is

%% fsolve for permeate concentration with linearized pore equation

global TempK Faraday R_id D_14_inf D_24 D_34
global VM_1 VM_2 VM_3 VM_4 
TempK = 25 + 273.15; % Temp in Kelvin

VM_1 = 1.8e-5;          % lysine molar volume [m3/mol]
VM_2 = 1.57e-5;         % Na molar volume [m3/mol]   * Calculated from rS
VM_3 = 4.47e-6;         % Cl molar volume [m3/mol]   * calcualted form rS
VM_4 = 1.8e-5;          % Water molar volume [m3/mol]
D_14_inf = 6.1e-11;     % D14 value under diluted conditions and IEP [m^2/s]
D_24 = 1.33E-9;         % D24 value under diluted conditions [m^2/s]
D_34 = 2.03E-9;         % D34 value under diluted conditions [m^2/s]
Faraday = 96485.336;       % Faraday constant [C/mol]
R_id = 8.314;            % R constant [J/(K mol)]

z_pol = 1.23e-5; 
r_AA = 3.48e-9; 
r_Na = 102e-12;   % ion size
r_Cl = 181e-12;   % ion size
rs = [r_AA; r_Na; r_Cl];
D = [ D_14_inf; D_24; D_34];
Vs = [VM_1; VM_2; VM_3];
rp = 9.75e-9; 
Le =  1.65e-4; % effective membrane thickness
[Kc, Kd, DP, phi] = HindranceFactor(rs, rp, D);
X_mem = -0.5;

C_ret_arb1 = [0.02 0.02 0.02];  % using average value of pH 7 0.15 M
C_ret_arb=[C_ret_arb1 1-sum(C_ret_arb1)];

Z_arb = [-1 1 -1];
Cp_guess_arb = [0.003 0.003 0.003]; % I changed the above line into the current line by removing C_ret

%% for arbitrary flux value
options = optimoptions('fsolve', 'Display','iter','algorithm', 'levenberg-marquardt','FunctionTolerance', 1e-4,'StepTolerance', 1e-4,'MaxFunctionEvaluation', 300);
Jv_arb = 5e-6;

[Cp_calc_arb, fval, exitflag] = fsolve(@(y)find_CP_MSNP( y, C_ret_arb,Jv_arb,Z_arb, z_pol,phi,X_mem), Cp_guess_arb(1:3),options); %from the line above

%have to find out why.
if exitflag < 1 
        iter = 0;
        mod = 0.01;
while xor(exitflag < 1, iter >50) 
        Cp_guess_arb = C_ret_arb .* [mod 1 1 1]; % can change all the value, this is just for demonstration purposes
        [Cp_calc_arb,fval, exitflag] = fsolve(@(y)find_CP_MSNP( y, C_ret_arb, Jv_arb, Z_arb, ...
                                z_pol), Cp_guess_arb(1:3), options);
        iter = iter +1
        mod = mod + 0.01;
end
end
Cp_calc_arb

Answer 1

[Matt J]

Why do I get different outputs for fsolve when I use different initial point?

Because you have local solutions. You need to choose your initial guess as close as possible to the global solution.

[Tadele]

Thank you very much!

Can I say that the approach I followed looks correct?

[Torsten]

Can I say that the approach I followed looks correct?

We don't know the problem you are trying to solve in a mathematical notation. We don't have your complete code. We don't know how your code works. So how could we tell ?

[Matt J]

Can I say that the approach I followed looks correct?

The use of global variables is discouraged. It would be better to use anonymous or nested functions as discussed in,

https://www.mathworks.com/help/optim/ug/passing-extra-parameters.html

[Tadele]

@Torsten and @Matt J

Thank you for your replies!

@Torsten It will really help me if you can look at the complete code I am trying to solve. The output of the model and experimental results should be similar. But there is a deviation between the two when I check. Therefore, your help will mean a lot to me.

The final aim of my code is to converge the difference between the calculated and guess value. But the calculated value should use the guess value first to get the calculated value. My codes are shown below one by one. I will start with the functions I called for the filnal solvers on number 6 and 7.

1. Function ''Hindrance factor''

function [Kc, Kd, DP, phi] = HindranceFactor(rs, rp, D_inf)

lambda = rs./rp;

phi=(1-lambda);

Kc = (2 - phi).*(1 + 0.054.*lambda - 0.988.*(lambda.^2) + 0.441.*(lambda.^3));

Kd = 1 - 2.3.*lambda + 1.154.*(lambda.^2) + 0.224.*(lambda.^3);

DP = Kd .* D_inf;

end

2. Differential equations to be solved using Ode15i

function dxdz = MS_AA_imp(z,x,xp, x_p, V, Z)

TempK = 25 + 273.15;    % Temp in Kelvin
Faraday = 96485.336;    % Faraday constant [C/mol]
R_id = 8.314;           % R constant [J/(K mol)]
VM_1 = 1.8e-5;          % AA molar volume of lysine [m3/mol] ---converted to m3--which was 18 cm3/mol
VM_2 = 1.57e-5;         % Na molar volume [m3/mol]   * Calculated from rS
VM_3 = 4.47e-6;         % Cl molar volume [m3/mol]   * calcualted form rS
VM_4 = 1.8e-5;          % Water molar volume [m3/mol]
D_24 = 1.33E-9;         % D24 value under diluted conditions [m^2/s]
D_34 = 2.03E-9;         % D34 value under diluted conditions [m^2/s]
h=7; 

hydration= x(1)/((x(1))*(1-x(1)*h)*(R_id*TempK));                 % new hydration formula
CT = x(1) / VM_1 + x(2)/ VM_2 + x(3) / VM_3 + (1-sum(x)) / VM_4;

D_14 =10e-9 *(0.665-0.393*x(1))/(1+x(1)* (h/(1-x(1)* h)));        % new D_14 formula
IStr = 0.5* (Z(2)^2 * x(2)/(VM_2 *1000) + Z(3)^2 * x(3)/(VM_3 *1000));
D_23 = ((D_24 + D_34)/2 * IStr ^(0.55) / (abs((Z(2)) * (Z(3)))^(2.3)));

J_1 = x_p(1) * V; 
J_2 = x_p(2) * V;
J_3 = x_p(3) * V;
J_4 = (1-sum(x_p)) * V;

%equations to be solved
dxdz = [hydration * xp(1) + x(1) * Z(1) * Faraday / (R_id .* TempK)*xp(5) - (J_4 * x(1) - J_1* x(4))/(D_14 * CT); 
 xp(2) + x(2) * Z(2) * Faraday*xp(5)  / ( R_id * TempK) - (J_3 * x(2) - J_2 * x(3))/ D_23 - (J_4 * x(2) - J_2 * x(4))/ (D_24 *CT);
 xp(3) + x(3) * Z(3) * Faraday*xp(5)  / ( R_id * TempK) - (J_2 * x(3) - J_3 * x(2))/ D_23 - (J_4 * x(3) - J_3 * x(4))/ (D_34* CT);   %D_32V replaced with D_    
 x(1) + x(2) + x(3) + x(4) - 1; 
 Z(1) * x(1) + Z(2)* x(2) + Z(3) * x(3)];
end

3. function ''partitionFirst". To be called (solved) by fsolve

function F = partitionFirst(C_entrance_MS,C_surf,Z,phi)
k=C_entrance_MS(1);       %concentration inside the membrane inlet
l=C_entrance_MS(2);
m=C_entrance_MS(3);
X_mem=C_entrance_MS(4);
VM_1 = 1.8e-5;          % lysine molar volume [m3/mol]
VM_2 = 1.57e-5;         % Na molar volume [m3/mol]   * Calculated from rS
VM_3 = 4.47e-6;         % Cl molar volume [m3/mol]   * calcualted form rS

F(1)= Z(1)*k + X_mem  + Z(2)*phi(2)*C_surf(2)*(k/(phi(1)*C_surf(1)))^(Z(2)/Z(1))...
    + Z(3)*phi(3)*C_surf(3)*(k/(phi(1)*C_surf(1)))^(Z(3)/Z(1));

F(2)= -l + phi(2)*C_surf(2)*(k/(phi(1)*C_surf(1)))^(Z(2)/Z(1));

F(3)= -m + phi(3)*C_surf(3)*(k/(phi(1)*C_surf(1)))^(Z(3)/Z(1));

F(4)= k*Z(1) + Z(2)*l + Z(3)*m + X_mem;
end

4. Function ''NernstPlanck_vol''. To be called (solved) by fsolve

function dydz = NernstPlanck_Vol(Co,C_entracne_MS, Z, Kc, DP) %linearized pore equations( derivated--can be found on other files)

xo_1 = Co(1); %concentration at the outlet of the membrane pore(end of membrane concentration)
xo_2 = Co(2);
xo_3 = Co(3);

TempK = 25 + 273.15;    % Temp in Kelvin
Faraday = 96485.336;    % Faraday constant [C/mol]
R_id = 8.314;           % R constant [J/(K mol)]
VM_1 = 1.8e-5;          % AA molar volume of lysine [m3/mol] ---converted to m3--which was 18 cm3/mol
VM_2 = 1.57e-5;         % Na molar volume [m3/mol]   * Calculated from rS
VM_3 = 4.47e-6;         % Cl molar volume [m3/mol]   * calcualted form rS
VM_4 = 1.8e-5;          % Water molar volume [m3/mol]
r_AA = 3.48e-9; 
r_Na = 102e-12;         % ion size
r_Cl = 181e-12;         % ion size
p_7 = 0.6e-3;     %Viscosity of pure water at 45C (Pa.s)(engineering toolbox.com)--denoted as p_7
deltaP = 5;       %pressure
D_im = 1.2e-9;    %diffusivitiy not yet worked out for all the solutes
Le =  1.65e-4;    % effective membrane thickness

C(1:3)=C_entracne_MS;
x_1 = C_entracne_MS(1); %correct by rreplaceing C() with C_entrnce()
x_2 = C_entracne_MS(2);
x_3 = C_entracne_MS(3);% here was s added by abel and I think it is an error

x_ave=(Co(1:3)+C(1:3))/2;

%corrected equtions
eq4 = ((x_ave(1)*Z(1)*r_AA^2*deltaP/(8*p_7*Le)*(Kc(1)-1)*R_id*TempK)/DP(1) - x_ave(1)*Z(1)*VM_1*deltaP/Le)/(Faraday*(Z(1)^2*x_ave(1)+Z(2)^2*x_ave(2)+ Z(3)^2*x_ave(3))) +...      
      ((x_ave(2)*Z(2)*r_Na^2*deltaP/(8*p_7*Le)*(Kc(2)-1)*R_id*TempK)/DP(2) - x_ave(2)*Z(2)*VM_2*deltaP/Le)/(Faraday*(Z(1)^2*x_ave(1)+Z(2)^2*x_ave(2)+ Z(3)^2*x_ave(3))) + ...
      ((x_ave(3)*Z(3)*r_Cl^2*deltaP/(8*p_7*Le)*(Kc(3)-1)*R_id*TempK)/DP(3) - x_ave(3)*Z(3)*VM_3*deltaP/Le)/(Faraday*(Z(1)^2*x_ave(1)+Z(2)^2*x_ave(2)+ Z(3)^2*x_ave(3)));

dydz(1) =x_1 - xo_1 + (x_ave(1)*(r_AA^2*deltaP/(8*p_7*Le)*(Kc(1)-1)*R_id*TempK)/DP(1) - x_ave(1)*VM_1*deltaP/Le - x_ave(1)*Z(1).*eq4)*Le; %assigning Elctric potential with number x(7)1

dydz(2) =x_2 - xo_2 + (x_ave(2)*(r_Na^2*deltaP/(8*p_7*Le)*(Kc(2)-1)*R_id*TempK)/DP(2) - x_ave(2)*VM_2*deltaP/Le - x_ave(2)*Z(2).*eq4)*Le; %assigning Elctric potential with number x(7)1

dydz(3) =x_3 - xo_3 + (x_ave(3)*(r_Cl^2*deltaP/(8*p_7*Le)*(Kc(3)-1)*R_id*TempK)/DP(3) - x_ave(3)*VM_2*deltaP/Le - x_ave(3)*Z(3).*eq4)*Le; %assigning Elctric potential with number x(7)1

end

5. function ''Partionsecond''. To be solved (called) by fsolve

function F=partitionsecond(Cp_calc,C_outlet,Z,phi)
kk=Cp_calc(1);
ll=Cp_calc(2);
mm=Cp_calc(3);

%bsed on electro neutrality in the permeate
F(1)= kk*Z(1) + Z(2)*C_outlet(2)/phi(2)*1/((C_outlet(1)/(phi(1)*kk))^(Z(2)/Z(1)))...
      + Z(3)*C_outlet(3)/phi(3)*1/((C_outlet(1)/(phi(1)*kk))^(Z(3)/Z(1))); %%% pationing equation for rearranging equating for Cp_cal
  
F(2)= -ll + C_outlet(2)/phi(2)*1/((C_outlet(1)/(phi(1)*kk))^(Z(2)/Z(1))); %partionig equation 2

F(3)= -mm + C_outlet(3)/phi(3)*1/((C_outlet(1)/(phi(1)*kk))^(Z(3)/Z(1))); %partionig equation 3
end

6. The above functions will be used to solve my problem. To do this, I used anonymous functions to use the output of one function as an input variable for the next function and so on. Below is the code for the cascaded functions.

function [sol] = find_CP_MSNP(Cp_guess, phi_start, Jv,Z, z_pol,phi,Kc,DP)
global VM_1 VM_2 VM_3 VM_4 
global TempK Faraday R_id D_14_inf D_24 D_34

%ode15i for MS_AA_imp
Cp_guess_MS = [Cp_guess 1- sum(Cp_guess)]; 
[y0, yp0]  = decic(@MS_AA_imp, 0, [phi_start 0], [0 0 0 0 0], [1 1 0 0 0], [],[],Cp_guess_MS, Jv, Z);
[zsol_MS, ysol_MS] = ode15i(@(z,x,xp)MS_AA_imp(z,x,xp ,Cp_guess_MS, Jv, Z), [0 z_pol], y0, yp0);

C_surf = ysol_MS(end,:)./[VM_1 VM_2 VM_3 VM_4 1];

%fsolve for partitionFirst
% define anonymous function 'f' for the partition function which accepts C_surf as an input parameter
zz=[0.00001 0.00001 0.00001 0.000001];
f = @(zz)partitionFirst(zz,C_surf(1:3),Z,phi);
C_entrance_pore = fsolve(f,zz);
C_entrance = C_entrance_pore.*[VM_1 VM_2 VM_3 1]; %to conver to mole fraction

%fsolve for NernstPlanck_Vol
% define anonymous function 'nernstPlanckAnonymous' for the NernstPlanck_Vol function which accepts C_entrance as a parameter
z=[0.00001 0.00001 0.00001];
nernstPlanckAnonymous = @(z)NernstPlanck_Vol(z,C_entrance(1:3),Z,Kc,DP);
C_end = fsolve(nernstPlanckAnonymous,z); %concentration on the end of membrane

% fsolve for partitionsecond uses the output of the above equation(C_end) as an
% input
secondpartitionAnonymous = @(z)partitionsecond(z, C_end(1:3),Z,phi); 
Cp_calc = fsolve(secondpartitionAnonymous,z);

sol = Cp_calc(1:3) - Cp_guess(1:3);

end

7. The last solver to solve the function on number 6 is shown below. Here the solver solves the function on number 6 (''find_CP_MSNP'') sothat the difference between calculated value (Cp_calc(1:3)) and guess value (Cp_guess(1:3)) to converge or approaches zero.

N.B. the guess value Cp_guess(1:3) is the starting point for the functions on number 6 (''find_CP_MSNP'') and used on the last point.

global TempK Faraday R_id D_14_inf D_24 D_34
global VM_1 VM_2 VM_3 VM_4 
TempK = 25 + 273.15; % Temp in Kelvin
VM_1 = 1.8e-5;          % lysine molar volume [m3/mol]
VM_2 = 1.57e-5;         % Na molar volume [m3/mol]   * Calculated from rS
VM_3 = 4.47e-6;         % Cl molar volume [m3/mol]   * calcualted form rS
VM_4 = 1.8e-5;          % Water molar volume [m3/mol]
D_14_inf = 6.1e-11;     % D14 value under diluted conditions and IEP [m^2/s]
D_24 = 1.33E-9;         % D24 value under diluted conditions [m^2/s]
D_34 = 2.03E-9;         % D34 value under diluted conditions [m^2/s]
Faraday = 96485.336;       % Faraday constant [C/mol]
R_id = 8.314;            % R constant [J/(K mol)]
z_pol = 1.23e-5; 
r_AA = 3.48e-9; 
r_Na = 102e-12;          % ion size
r_Cl = 181e-12;          % ion size
rs = [r_AA r_Na r_Cl];
D = [ D_14_inf D_24 D_34];
Vs = [VM_1 VM_2 VM_3];
rp = 9.75e-9; 
Le =  1.65e-4;           % effective membrane thickness
[Kc, Kd, DP, phi] = HindranceFactor(rs, rp, D)

C_ret_arb1 = [0.02 0.02 0.02];  % using average value of pH 7 0.15 M
C_ret_arb=[C_ret_arb1 1-sum(C_ret_arb1)];

Z_arb = [-1 1 -1];
Cp_guess_arb = [0.025  0.033 0.023]; % I changed the above line into the current line by removing C_ret

%% for arbitrary flux value
options = optimoptions('fsolve', 'Display','iter','algorithm', 'levenberg-marquardt','FunctionTolerance', 1e-4,'StepTolerance', 1e-4,'MaxFunctionEvaluation', 300);
Jv_arb = 5e-6;

[Cp_calc_arb, fval, exitflag] = fsolve(@(y)find_CP_MSNP( y, C_ret_arb,Jv_arb,Z_arb, z_pol,phi,Kc,DP), Cp_guess_arb(1:3),options); %from the line above

%have to find out why.
if exitflag < 1 
        iter = 0;
        mod = 0.01;
while xor(exitflag < 1, iter >50) 
        Cp_guess_arb = C_ret_arb .* [mod 1 1 1]; % can change all the value, this is just for demonstration purposes
        [Cp_calc_arb,fval, exitflag] = fsolve(@(y)find_CP_MSNP( y, C_ret_arb, Jv_arb, Z_arb, ...
                                z_pol,phi,Kc,DP), Cp_guess_arb(1:3), options);
        iter = iter +1
        mod = mod + 0.01;
end
end
Cp_calc_arb

I hope the codes are clear and you will be able to help me.

Thank you very much in advance.