why fsolve cant find the solutoin?

Question

Mohiedin Bagheri 2023-1-17

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1895300-why-fsolve-cant-find-the-solutoin

评论： Star Strider 2023-1-17

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2023-1-17

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1895300-why-fsolve-cant-find-the-solutoin#answer_1150205

编辑：Star Strider 2023-1-17

在 MATLAB Online 中打开

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

.