This question suffers from several misunderstandings about optimization tools, several common problems.
First, think of an optimizer as a blindman, placed somewhere on the face of the earth. Give him only his cane, then task him with finding the lowest spot on earth. Yes, I hope you will give him scuba gear too, as that spot will be deep in a Pacific ocean trench. But the point is, if you set him down to start with in the Himalayas, do you exeect him to find the Marianas trench in the Pacific? As bad, supose you set him down near the Dead Sea? He will find a low spot, but then decide that no direction leads him down farther, so he will stop searching.
An optimization tool is the same. If you start it in the wrong place, then what do you expect? It looks aroun, then may get lost if you give it different starting values. All is well and good when a problem has only ONE local minimizer, so that solution is also the globally best solution. But many problems do not have this property. So change the starting values, and what do you expect? Good starting values often matter.
Next, problems that vary over many powers of 10 in the paameters are nasty things. Worse, your problem goes to hell if the optimizer looks even slightly below zero, yet your parameters are very near zero in absolute terms.
The good thing is this is easily repaired. Work in terms of log to the base 10. So define the unknowns as log10 of your unknowns.
Your initial guess now becomes
r0 = [-5 -6 -7 -9] - 1;
% Initial guess
Inside your objective, the very first thing you will do is raise 10 to the power of r.
function F = myfun(r,pH,p)
R = 10.^r;
c_CO2 = R(1);
c_H2CO3 = R(2);
c_HCO3 = R(3);
c_CO3 = R(4);
...
But as far as the optimizer is concerned, now all the unknowns have the same order of magnitude, and things are now well scaled. As well, it never has to worry about numbers going less than zero, since the transformation handles everything.
Oh, and trying to set TolFun to 1e-20? Silly. You can't achieve that much accuracy in a problem like this.
So, first, I'll hack your objective function a bit, just to avoid needing to pass in those parameters. I am feeling lazy today.
function F = myfun(r)
p.k_H = 29.76; % Henry constant
p.K_a1 = 2.5*10^-4; % Equilibrium constant 1 (H2CO3 <-> HCO3- + H+)
p.K_a2 = 4.69*10^-11; % Equilibrium constant 2 (HCO3- <-> CO3-- + H+)
p.K_h = 10^-2.8; % Hydration constant (H2CO3/CO2)
pH = 4.92; % pH-value
R = 10.^r;
c_CO2 = R(1);
c_H2CO3 = R(2);
c_HCO3 = R(3);
c_CO3 = R(4);
c_H = 10^-pH;
F(1) = 10^-14/c_H + c_HCO3 + 2*c_CO3 - c_H; % Charge neutrality
F(2) = c_H2CO3/c_CO2 - p.K_h; % Hydration
F(3) = c_H*c_CO3/c_HCO3 - p.K_a2; % Eq2
F(4) = c_H*c_HCO3/c_H2CO3 - p.K_a1; % Eq1
Next, test it!!!!!!!! Verify that it does something reasonable.
myfun([-6 -7 -8 -9])
ans =
-1.201e-05 0.098415 1.2022e-06 -0.0002488
So it produces numbers, 4 of them as expected, varying over 4 orders of magnitude. That in itself will be highly problematic. Why? A tool like lsqnonlin tries to minimize the sum of squares of the numbers it sees. So if one of those values is 4 powers of 10 larger than the others, it can essentially ignore the other three. All that matters is the biggest one.
So merely using lsqnonlin here MAY be a highly suspect thing to do. And settign TolFun at 1e-20? I'd bet it is way too optimistic. But for that tight a tolerance, you will need to increase the maximum function evals allowed. (And learn to use scientific notation with your numbers.)
Next, set the display parameter to iter, so we can see if it is converging. ALWAYS do this the first time you run a specific optimization.
options = optimset('TolFun',1e-15,'MaxFunEvals',5000,'MaxIter',1000,'Display','iter');
We can now try lsqnonlin, but there is no longer a need for a lower bound, because we are working in a log space instead.
r0 = [-6 -7 -8 -9];
val = lsqnonlin(@myfun,r0,[],[],options)
Norm of First-order
Iteration Func-count f(x) step optimality
0 5 0.0096856 0.0227
1 10 0.00129798 10 0.00308
2 15 0.00129798 5.31272 0.00308
3 20 0.000173209 1.32818 0.000421
4 25 2.00996e-05 0.948856 5.85e-05
5 30 1.78119e-06 1.02082 8.18e-06
6 35 8.57658e-08 1.05491 1.01e-06
7 40 3.71167e-09 1.06807 5.56e-08
8 45 3.59717e-10 1.00663 1.35e-09
9 50 3.54168e-10 0.00409291 6.69e-10
10 55 3.48651e-10 0.00591474 1.33e-09
11 60 3.43293e-10 0.00399236 6.4e-10
12 65 3.37891e-10 0.00597722 1.31e-09
13 70 3.32645e-10 0.00389634 6.26e-10
14 75 2.0284e-11 0.828326 1.77e-10
15 80 2.01912e-11 0.000524991 1.06e-10
16 85 2.00988e-11 0.000775195 1.76e-10
17 90 2.00068e-11 0.00052201 1.05e-10
18 95 1.99152e-11 0.000772959 1.75e-10
19 100 1.98242e-11 0.000519053 1.04e-10
20 105 1.97335e-11 0.000770741 1.75e-10
21 110 1.96433e-11 0.000516122 1.04e-10
22 115 1.95534e-11 0.000768542 1.74e-10
...
980 4905 1.12224e-14 8.67351e-05 7.73e-12
981 4910 1.10723e-14 1.77985e-05 1.27e-12
982 4915 1.09241e-14 8.56892e-05 7.63e-12
983 4920 1.07778e-14 1.75701e-05 1.25e-12
984 4925 1.06333e-14 8.46533e-05 7.53e-12
985 4930 1.04908e-14 1.73441e-05 1.24e-12
986 4935 1.035e-14 8.36273e-05 7.43e-12
987 4940 1.02111e-14 1.71207e-05 1.22e-12
988 4945 1.00739e-14 8.26112e-05 7.34e-12
989 4950 9.93853e-15 1.68998e-05 1.2e-12
990 4955 9.80485e-15 8.1605e-05 7.24e-12
991 4960 9.67296e-15 1.66814e-05 1.19e-12
992 4965 9.54272e-15 8.06085e-05 7.15e-12
993 4970 9.41422e-15 1.64655e-05 1.17e-12
994 4975 9.28731e-15 7.96218e-05 7.06e-12
995 4980 9.16212e-15 1.6252e-05 1.15e-12
996 4985 9.03848e-15 7.86448e-05 6.96e-12
997 4990 8.91651e-15 1.60409e-05 1.14e-12
998 4995 8.79606e-15 7.76776e-05 6.87e-12
999 5000 8.67723e-15 1.58322e-05 1.12e-12
1000 5005 8.55989e-15 7.672e-05 6.78e-12
val =
-3.4349 -6.2349 -4.9169 -8.481
I skipped the middle 900 or so iterations in that list.
What did that get us to?
myfun(val)
ans =
2.9286e-06 2.623e-09 2.6118e-09 4.1699e-08
Trying to set the tolerance even that low is surely twiddling with how many angels can fit on the head of a pin. But lsqnonlin seems to work reasonably here. The biggest improvement was surely moving to the logs of your parameters, as that made the problem well scaled in the eyes of lsqnonlin.
