why fsolve cant find the solutoin?

1 view (last 30 days)
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

Accepted Answer

Star Strider
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
.
  6 Comments
Mohiedin Bagheri
Mohiedin Bagheri on 17 Jan 2023
Thank you @Star Strider, I aprpeciate it
All the best

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!