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

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.

Answer 1

[Alan Weiss]

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

Answer 2

[Mike Croucher]

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?