why fsolve cant find the solutoin?

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

Answer 1

[Star Strider]

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

.

[Mohiedin Bagheri]

Thank you so much @Star Strider. I really appreciate your guidance.

[Star Strider]

As always, my pleasure!

[Mohiedin Bagheri]

Dear Star Strider,

This question that you gently answered was built based on a novel approach I developed in my PhD dissertation The following lines are very confidential for me, would you be able to modify or change the following lines so it be hidden in the answer you provided?

I didn't notice I shared the whole code including these lines when I shared the .m file in public here yesterday.

I really appreciate your understanding.

I appreciate if you could hide these sections below or modify it someway in your answer:

%Potential-independent kinetic factors obtained from transient

%Potential-dependent kinetic factors

%% Electrolyte (Solution Thermodynamics):

Thanks a lot Star Strider again

[Star Strider]

I understand.

I deleted the requested lilnes. I kept a copy of the entire code in the event I need to refer to it.

[Mohiedin Bagheri]

Thank you @Star Strider, I aprpeciate it

All the best

[Star Strider]

As always, my pleasure!