why fsolve cant find the solutoin?
1 view (last 30 days)
Show older comments
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
0 Comments
Accepted Answer
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
Ex = interp1(ic, E, icx) % 'E' Correspnding Coordinate
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);
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
More Answers (0)
See Also
Categories
Find more on Startup and Shutdown in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!