fsolve gives the same value for different parameters of the nonlinear equation

3 views (last 30 days)
Fabian Urias
Fabian Urias on 29 Nov 2021
Answered: Mike Croucher on 30 Nov 2021
fsolve gives me the same value for H, changing the parameters in the function
function f=funH(x,pmt)
H=x(1);
% Parmetrs
C_TIC=pmt(1);
c_A=pmt(2);
T=pmt(3);
KC=exp(31.78-8448/T)/10^3;
K1=10^-(-14.8435+0.032786*T+3404.71/T);
K2=10^-(-6.4980+0.02379*T+2902.39/T);
KW=exp(-(148.9802-13847.26/T+23.6521*log(T)));
HCO3=K1*KC*H*C_TIC/(H^2+KC*H^2+KC*K1*H+KC*K1*K2);
CO3=KC*K1*K2*C_TIC/(H^2+KC*H^2+KC*K1*H+KC*K1*K2);
OH=KW/H;
f=OH+HCO3+2*CO3-H-c_A;
The only parameter I change is C_TIC, with that it is supposed to give me different values of H when changing it, but I always get the same value of -100.
The main code is the following:
cc=[0
0.0273130929895407
0.0544594523216784
0.0814352072934854
0.108239495044338
0.134870346309588
0.161328561502616
0.187615205950685
0.213730223059828
0.239675750650148
0.265453667872101
0.291072309701383
0.316550740405870
0.341867505185996
0.367022604041762
0.392016036973168
0.416847803980213
0.441517905062897
0.466026340221222
0.490373109455186
0.514569485192549
0.538660214173914
0.562602127209587
0.586395224299569
0.610039505443860
0.633534970642459
0.656881619895367
0.680079453202584
0.703128470564109
0.726034917926889
0.748861064536146
0.771552936694940
0.794110534403270
0.816533857661137
0.838822906468540
0.860977680825479
0.882998180731955
0.904884406187968
0.926636357193516
0.948325704755969
0.969895136668155
0.991344523936537
1.01267386656111
1.03388316454189
1.05497241787886
1.07594162657202
1.09679079062138
1.11751991002693
1.13816455304136
1.15870592446615
1.17914222665592
1.19947488725981
1.21970533392692
1.23983499430640
1.25986529604738
1.27979766679897
1.29963353421031
1.31937432593052
1.33902146960873
1.35857639289408
1.37804052343568
1.39741528888267
1.41669806121460
1.43588905353830
1.45499320028932
1.47401181822714
1.49294622411125
1.51179773470111
1.53056766675622
1.54925733703604
1.56786806230005
1.58640115930774
1.60485794481859
1.62323973559206
1.64154784838764
1.65978359996482
1.67794696939914
1.69602737753052
1.71403647939947
1.73197540500091
1.74984528432973
1.76764724738085
1.78538242414919
1.80305194462963
1.82065693881711
1.83819853670652
1.85567786829278
1.87309606357079
1.89045425253547
1.90775356518172
1.92499513150445
1.94216926722203
1.95928065498771
1.97633434852469
1.99333127053802
2.01027234373278
2.02715849081404
2.04399063448686
2.06076969745630
2.07749660242745
2.09417227210535
2.11079762919509
2.12737359640173
2.14390109643032
2.16038105198595
2.17681048812212
2.19318348188068
2.20950886241089
2.22578737175450
2.24201975195323
2.25820674504881
2.27434909308298
2.29044753809748
2.30650282213403
2.32251568723436
2.33848687544021
2.35441712879332
2.37030718933540
2.38615779910821
2.40196970015346
2.41773631483513
2.43346296730304
2.44915087405706
2.46480054407189
2.48041247661213
2.49598716123227
2.51152507777669
2.52702669637969
2.54249247746544
2.55792287174802
2.57331832023140
2.58867925420945
2.60400609526592
2.61929925527448
2.63455942529360
2.64978714797925
2.66498259641862
2.68014617598014
2.69527828414320
2.71037931049813
2.72544963674621
2.74048963669968
2.75549967628172
2.77048011352645
2.78543129857897];
options = optimset('Display','off');
for i=1:length(cc)
[xx]=fsolve(@(x)funH (x,[cc(i) 100 303.15]),[5], options);
H(i)=xx(1);
end
I'm sorry for including so many parameters, I think it would be easier for you guys to add it so than include it in an excel file.

Answers (2)

Alan Weiss
Alan Weiss on 30 Nov 2021
Your function funH returns the following iinternal values for your initial point x0 = 5 and the other parameters pmt = [1.0273130929895407,100,303.15]:
f=OH+HCO3+2*CO3-H-c_A where
OH = 5.5e-105
HCO3 = 4.6e-9
CO3 = 4.75e-20
H = x = 5
c_A = 100
In other words, all of the portions of your output f are tiny except for x and c_A = 100. So the solution is always going to be very near x = -100, which is where -x -c_A = 0.
Alan Weiss
MATLAB mathematical toolbox documentation

Mike Croucher
Mike Croucher on 30 Nov 2021
It seems that chaging the paramet C_TIC from 0 to 2.78 doesn't affect the outpur much. I changed your function slightly so that its vectorised. This makes it easier to plot
function f=funH(x,pmt)
H=x;
% Parmetrs
C_TIC=pmt(1);
c_A=pmt(2);
T=pmt(3);
KC=exp(31.78-8448/T)/10^3;
K1=10^-(-14.8435+0.032786*T+3404.71/T);
K2=10^-(-6.4980+0.02379*T+2902.39/T);
KW=exp(-(148.9802-13847.26/T+23.6521*log(T)));
HCO3=K1.*KC.*H.*C_TIC./(H.^2+KC*H.^2+KC.*K1.*H+KC*K1*K2);
CO3=KC*K1*K2*C_TIC./(H.^2+KC*H.^2+KC.*K1.*H+KC.*K1.*K2);
OH=KW./H;
f=OH+HCO3+2*CO3-H-c_A;
end
Now let's create two versions of this, one with a C_TIC of 0 and one with a C_TIC of 2.78 which represent the extremes of your cc list
func1 = @(x)funH (x,[0 100 303.15]);
func2 = @(x)funH (x,[2,78 100 303.15]);
we can now plot these. Here's func1
x = linspace(-110,-90,50);
y1 = func1(x); % recall this has a C_TIC of 0
plot(x,y1)
and next we we have func2
y2 = func2(x);
plot(x,y1);
These look identical to my eyes! Made me wonder if there was any difference at all. Let's plot the difference and see
plot(x,y1-y2);
So they are different but look at the Y axis of the difference plot. These differences are of the order of 10^-10 when the absolute values across the same x-range are of the order of go from -10 to 10.
When you look at your output array, H, which contains all the results of the fsolve command it does indeed look like they are all the same
H =
Columns 1 through 12
-100.0000 -100.0000 -100.0000 -100.0000 -100.0000
but if you plot it you'll see that they are different
plot(H)
In short, your code works fine but maybe you might want to rescale your function somehow?

Tags

Products


Release

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by