Trying to use fsolve works over lower part of range (values < 10) but fails in upper part of range (values > 10)

Running the following partial script gives me results as expected:

Kw = 1*10^-14;
Ca = 0.1;
Va = 10;
pKa = 4.757;
Ka = 10^-pKa;
Cb = 0.1;
pKb = 0.2;
Kb = 10^-pKb;

pHv = (2.9:0.2:12.5);
Vbv = zeros(1,numel(pHv));
for k = 1:numel(pHv)
    H = 10^-pHv(k);
    Vbv(k) = Va * ((Ca/(1+H/Ka) - H + Kw/H)/(Cb/(1+Kw/(H*Kb )) + H - Kw/H ));
end
plot(Vbv,pHv)

Sample results

pH = 2.9 Vb = 0.011058027

pH = 8.1 Vb = 9.995732731

pH = 9.3 Vb = 10.00402116

pH = 12.5 Vb = 20.69392396

All of the above is known to be correct. In practice, I want to perform the calculation in reverse. In other words calculate pH for given values of Vb.

The following has been suggested (same initial values as above) but I don't really understand it:

Vbfcn = @(Ca,Cb,Ka,Kb,Kw,Va,H) Va.*((Ca./(1+H./Ka) - H + Kw./H) ./ (Cb./(1+Kw./(H*Kb)) + H - Kw./H));

In the line above, I don't understand the purpose of the full stops (e.g. Cb.) or the @

Is @(Ca,Cb,Ka,Kb,Kw,Va,H) a way to specify a list of parameters? and is Cb. (for example) a sort of placeholder for the Cb parameter?

Vbv = (0:0.2:20);

for m = 1:numel(Vbv)

Hv(m) = fsolve(@(H)norm(Vbw(m) - Vbfcn(Ca,Cb,Ka,Kb,Kw,Va,H)), 1E-11);

end

I am not clear what is going on in the fsolve line but I am guessing that 1E-11 is possibly an initial guess value? What is norm doing?

This seems to work for values of Vb < 10 but not for values of Vb > 10 so I am wondering if a different initial guess is needed for those values

I have tried the following while testing:

H = fsolve(@(H)norm(0.011058027 - Vbfcn(Ca,Cb,Ka,Kb,Kw,Va,H)), 1E-11)

pH = -log10(H)

Result H = 0.0013, pH = 2.9000 CORRECT

H = fsolve(@(H)norm(9.995732731 - Vbfcn(Ca,Cb,Ka,Kb,Kw,Va,H)), 1E-11)

pH = -log10(H)

Result H = 7.7682e-09, pH = 8.1097 NEAR ENOUGH CORRECT?

But here things go awry:

H = fsolve(@(H)norm(10.00402116 - Vbfcn(Ca,Cb,Ka,Kb,Kw,Va,H)), 1E-11)

pH = -log10(H)

Result H = -6.6663e-09, pH = 8.1761 - 1.3644i (Should be 9.3)

Answer 1

[Torsten]

It's very important that you have a good initial guess for pHv given Vbv.

So looping over the complete range of Vbv and taking the solution of the last step as initial guess for the next step is a good idea.

Vbv = 0:0.2:20;
pHv0 = 2;
options = optimset('Display','off','TolFun',1e-12,'TolX',1e-12);
for i = 1:numel(Vbv)
    fun_actual = @(x)fun(x) - Vbv(i);
    pHv(i) = fsolve(fun_actual,pHv0,options);
    error(i) = abs(fun_actual(pHv(i)));
    pHv0 = pHv(i);
end
figure(1)
plot(Vbv,pHv)

figure(2)
plot(Vbv,error)

function Vbv = fun(pHv)
  Kw = 1*10^-14;
  Ca = 0.1;
  Va = 10;
  pKa = 4.757;
  Ka = 10^-pKa;
  Cb = 0.1;
  pKb = 0.2;
  Kb = 10^-pKb;
  H = 10^-pHv;
  Vbv = Va * ((Ca/(1+H/Ka) - H + Kw/H)/(Cb/(1+Kw/(H*Kb )) + H - Kw/H ));
end

Answer 2

[Walter Roberson]

 Vbfcn = @(Ca,Cb,Ka,Kb,Kw,Va,H) Va.*((Ca./(1+H./Ka) - H + Kw./H) ./ (Cb./(1+Kw./(H*Kb)) + H - Kw./H));

In the line above, I don't understand the purpose of the full stops (e.g. Cb.) or the @

Is @(Ca,Cb,Ka,Kb,Kw,Va,H) a way to specify a list of parameters?

Yes, @(VARIABLES) specifies nicknames to use for whatever value is passed in to the corresponding position. That anonymous function definition is nearly exactly the same as

 Vbfcn = @(varargin) varargin{6}.*((varargin{1}./(1+varargin{7}./varargin{3}) - varargin{7} + varargin{5}./varargin{7}) ./ (varargin{2}./(1+varargin{5}./(varargin{7}*varargin{4})) + varargin{7} - varargin{5}./varargin{7}));

that is, the function receives a parameter list, and substitutes the value passed in for each position in place of the corresponding name. The names used are not important to MATLAB: it is all just a matching process, name inside the body gets compiled in as a reference to a positional argument.

The difference between the two different potential definitions are:

• using names is easier to read
• using names is easier to debug
• in the uncommon case where a function is called that uses inputname then varargin{NUMBER} would come out as '' (empty character vector) but as a character string representing the name of the anonymous function parameter in the case names are used
• the varargin version does not test for too many input parameters passed in, but the version with named parameters does

@ can be followed by:

• an unqualified unindexed variable name, such as @sin . In that case MATLAB looks through its resolution priority list to find a function with that name, and creates a reference to that function. It is completely legal in MATLAB to have a local or private function with the same name as a function used elsewhere; the @ handle created will refer specifically to the function with that name that you would get if you were to use that function name at that point, and that reference will hold even if the function handle is returned out of context to a place that cannot normally invoke the function that is being referenced. This creates a "simple" function handle
• a qualified name of a function in a package. For example @NOAA.atmosphere.scatter to reference +NOAA/+atmosphere/scatter.m . Same rules as above, including creating a "simple" function handle
• the name of a class
• A list of unindexed names inside () . In this case it creates an anonymous function, with parameter names as described above

and is Cb. (for example) a sort of placeholder for the Cb parameter?

MATLAB ran out of single-character operator names . It needs operators to indicate algebraic matrix multiplication ("inner product"), and matrix power (A^2 meaning A*A where * is inner product), and left and right fitting operations similar to inv(A)*B and A*inv(B); and it also needs element-by-element multiplication, and element-by-element power, and element-by-element right division (element by element left division is rare). MATLAB chose to use single character operators for the matrix versions (like inner product) and to use period before the similar operation for element-by-element work. So H/Ka would mean approximately inner_product(H, pinv(Ka)) and H./Ka would mean element_by_element_division(H, Ka)

[Walter Roberson]

     Hv(m) = fsolve(@(H)norm(Vbw(m) - Vbfcn(Ca,Cb,Ka,Kb,Kw,Va,H)), 1E-11);

I am not clear what is going on in the fsolve line but I am guessing that 1E-11 is possibly an initial guess value?

Yes, exactly.

What is norm doing?

When you are trying to balance two expressions, f(x) == g(x) when you can take f(x)-g(x) and aim for 0. In the case where both expressions are scalar and certain to be real-valued then you can just leave it like that and fsolve() can deal with it. In the case where either side might be vector valued then f(x)-g(x) would be vector valued; it happens that fsolve() can deal with that case too -- but not all solvers can handle that, so you tend to instead look for the sum of squares being 0. For real values, sqrt() of the sum of squares is norm()

In the case where either side might potentially have complex values then fsolve() becomes less reliable, and you might need to give it the hint of using norm() .

In this situation, norm() is probably not needed for fsolve(), but some people prefer to express the target mathematically so that they do not have to concentrate on the details of exactly what the solver can or cannot handle.

[Paul]

Well Walter, I asked for some explanantion and you delivered in spades! The @ explanation I can grasp fairly easily but I will need to spend a bit more time to get my head around the rest of it. For example I would probably like to explore the effects of omitting the periods but I'm sure I get the gist of it.

Many thanks, much appreciated.

[Walter Roberson]

M = magic(3)

M = 3×3

8 1 6 3 5 7 4 9 2

M * M  %inner product, algebraic matrix multiplication

ans = 3×3

91 67 67 67 91 67 67 67 91

M .* M  %element by element multiplication

ans = 3×3

64 1 36 9 25 49 16 81 4

M^2    %matrix power, same as M*M (inner product)

ans = 3×3

91 67 67 67 91 67 67 67 91

M.^2   %element by element squaring

ans = 3×3

64 1 36 9 25 49 16 81 4

M/M   %matrix right division, like M * pinv(M)

ans = 3×3

1 0 0 0 1 0 0 0 1

M./M  %element by element division

ans = 3×3

1 1 1 1 1 1 1 1 1

N = [10;20;30]

N = 3×1

10 20 30

M * N          %inner product, matrix multiplication

ans = 3×1

280 340 280

M .* N         %element-by-element with implicit expansion

ans = 3×3

80 10 60 60 100 140 120 270 60

M ./ N         %element-by-element with implicit expansion

ans = 3×3

0.8000 0.1000 0.6000 0.1500 0.2500 0.3500 0.1333 0.3000 0.0667

M / N          %matrix right division

Error using /
Matrix dimensions must agree.

[Walter Roberson]

When either side of A*B is scalar the result is the same as A.*B -- the scalar is multiplied by each element in the other matrix

When the right side of A/B is scalar the result is the same as A./B -- the scalar is divided into each element of A.

When the left side of A/B is scalar but the right side is not, then it is an error if B does not have exactly one column. If B does have exactly one column and A is scalar then A/B becomes a fitting operation

format long g
A = 3; B = [2; 5; 9];
A/B

ans = 1×3

0 0 0.333333333333333

When A and B are scalars then A^B is the same as A.^B . But you can see that A^B is a completely different operation when B is non-scalar

format long g
A = 3

A = 
     3

A^[1 2; 3 4]

ans = 2×2

87.8863782333136 127.119793533935 190.679690300902 278.566068534216

A.^[1 2;3 4]

ans = 2×2

3 9 27 81

Generally speaking, you should use .* and .^ and ./ for multiplication and power and division, unless you specifically mean matrix operations, with the possible exception of matrix * constant or constant * matrix or matrix / constant as an operation such as M/2 is pretty clear .

[Walter Roberson]

in MATLAB, the least-squared solution to A*x = b with known A and B is expressed as A\b and the solution to x*A = b is expressed as b/A . A\b is mathematically similar to pinv(A)*b and b/A is mathematically similar to b*pinv(A)

But b/A is not element-by-element division of the elements of a and the elements of A.

A common mistake people make is something like

x = 0:10;
x / (x^2 + 1)

Error using ^
Incorrect dimensions for raising a matrix to a power. Check that the matrix is square and the power is a scalar. To operate on each element of the matrix individually, use POWER (.^) for elementwise
power.

This has two problems: x^2 is x*x which is inner-product (algebraic matrix multiplication), which requires that the number of columns on the left side of x*x be equal to the number of rows on the right side of x*x . When x is a vector, that can never be true. x has to be a square matrix (including scalar) for x^2 to be defined at all, and then it is inner product of a matrix with itself rather than being element-by-element squaring. So the x^2 will not work for vector x.

If you fix up to x / (x.^2 + 1)

then you have row vector x matrix divide by row vector (x.^2 + 1) and that is an error.

In such cases people almost always want x ./ (x.^2 + 1) so that element-by-element operations are done.

If you want each element of a matrix to be treated individually, use the dot form of the operator ... which is most of the time for most code. You might get away with using the un-dotted operators when you use scalars, but it is a bad habit (and one that makes it more difficult to extend the code to handle non-scalars.)

Answer 3

[Paul]

Torsten, I am immensely grateful BUT I haven't actually learnt very much. If some of the questions in my original post could be answered I would at least know if I was thinking along the right lines. Also, if you could please explain or walk me through your code I would be better placed to understand how and why it works. I could do with a better understanding of the following lines:

options = optimset('Display','off','TolFun',1e-12,'TolX',1e-12);
fun_actual = @(x)fun(x) - Vbv(i);
pHv(i) = fsolve(fun_actual,pHv0,options);

Lastly, but perhaps less importantly, I am not clear what the error plot is supposed to tell me ...

[Torsten]

I think it would be more helpful for you if you use the possibility to visit the MATLAB Onramp course:

https://de.mathworks.com/learn/tutorials/matlab-onramp.html

Also reading the "fsolve" documentation will answer many of your questions.

The error plot is a plot of the values

abs(Vbv - (Va * ((Ca/(1+H/Ka) - H + Kw/H)/(Cb/(1+Kw/(H*Kb )) + H - Kw/H ))))

and shows how good the determined value for H for given Vbv satisfies

Vbv - (Va * ((Ca/(1+H/Ka) - H + Kw/H)/(Cb/(1+Kw/(H*Kb )) + H - Kw/H ))) = 0

[Paul]

Thank you Torsten, the Onramp course was extremely helpful (if I had known about it sooner I wouldn't have asked so many dumb questions). I also went over the fsolve documentation again, and it now makes more sense. The documentation on Tolerances and Stopping Criteria also helped me make sense of the optimset function.

[Torsten]

Thank you Torsten, the Onramp course was extremely helpful (if I had known about it sooner I wouldn't have asked so many dumb questions). I also went over the fsolve documentation again, and it now makes more sense. The documentation on Tolerances and Stopping Criteria also helped me make sense of the optimset function.

Nice to hear. Furthermore good luck.