# why fsolve cant find the solutoin?

3 views (last 30 days)
Mohiedin Bagheri on 17 Jan 2023
Commented: Star Strider on 17 Jan 2023
Hello,
I am trying to use fsolve to find intersection of two functions (intersection of two functions 'ia' and 'ic') . please see the attached code.
I ploted two functions (ia, and ic) and I can see the estimate where they cross, but the solution that fsolve gives at the end is totally wrong. It doesn't work correct! Would you please guide me what is wrong with my code and help me fixate the wrong lines of the code?
Thank you very much in advance

Star Strider on 17 Jan 2023
Edited: Star Strider on 17 Jan 2023
I cannot understand ther fsolve function (or the call to it), and in any event, it is not necessary since it is possible calculate the x-coordinate of the ‘ia-ic’ intersection with a simple interp1 call, as I did here. The appropriate value for ‘E’ can be calculated similarly —
%% 1- insert the operationl conditions:
pH = 4 ; T = 298 ; F = 96485;
E = linspace(-0.8,0,600);
%% Electrode (Surface Kinetic):
% — SEVERAL LINES IN THE CODE ARE DELETED FOR CONFIDENTIALITY REASONS, AS REQUESTED BY Mohiedin Bagheri —
ic = 1e2*10.^(E);
%i_tot = abs(ia-ic);
icx = interp1(ia-ic,ic,0) % 'ic' Intersection Coordinate
icx = 40.7639
Ex = interp1(ic, E, icx) % 'E' Correspnding Coordinate
Ex = -0.3897
figure(1)
plot(ia,E)
hold on
plot(ic,E)
%plot(i_tot,E)
set(gca, 'XScale', 'log')
set(gca, 'YLim', [-1.5 0])
set(gca, 'XLim', [0.001 10000])
ylabel('E vs. SHE (V)')
xlabel('Current density (A/m2)')
legend('ia','ic', 'Location','best')
xl = xline(icx, '-m', 'ic ia Intersection', 'DisplayName',sprintf('ic = ia = %.2f',icx));
xl.LabelVerticalAlignment = 'middle';
xl.LabelHorizontalAlignment = 'center';
yl = yline(Ex, '-m', 'E Intersection', 'DisplayName',sprintf('E = %.4f',Ex));
yl.LabelVerticalAlignment = 'middle';
yl.LabelHorizontalAlignment = 'center'; iEo=[0.01 -0.5];
options = optimset('Largescale', 'off','Display','off');
iE=fsolve(@Corr,iEo,options,k01,k02,k03,k04,b1,b2,b3,b4,k0_3a,k0_3t,b_3a,b_3t,F);
disp(iE);
-0.0000 -8.3083
function Z=Corr(iE,k01,k02,k03,k04,b1,b2,b3,b4,k0_3a,k0_3t,b_3a,b_3t,F);k1 = k01*exp((2.3/b1).*iE(2)); k2 = k02*exp((2.3/b2).*iE(2)); k3 = k03*exp((2.3/b3).*iE(2));
k4 = k04*exp((2.3/b4).*iE(2)); k_3 = k0_3a*exp((2.3/b_3a).*iE(2))+k0_3t*exp((2.3/b_3t).*iE(2));
theta_stst_1 = (k1.*k_3)./(k1.*k3+k1.*k_3+k2.*k_3); i1 = 2*F*k2.*theta_stst_1;
theta_stst_2 = (k1.*k3)./(k1.*k3+k1.*k_3+k2.*k_3); i2 = 2*F*k4.*theta_stst_2; ia = i1+i2;
%% Electrolyte (Solution Thermodynamics):
Z=zeros(2,1);
Z(1)= iE(1)-(i1+i2);
Z(2)= iE(1)-1e2*10.^(iE(2));
end
.
Star Strider on 17 Jan 2023
As always, my pleasure!