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

Question

Paul 2022-11-25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1863318-trying-to-use-fsolve-works-over-lower-part-of-range-values-10-but-fails-in-upper-part-of-range

评论： Torsten 2022-11-27

采纳的回答： Torsten

在 MATLAB Online 中打开

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)

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Torsten 2022-11-25

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1863318-trying-to-use-fsolve-works-over-lower-part-of-range-values-10-but-fails-in-upper-part-of-range#answer_1112023

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

Walter Roberson 2022-11-25

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1863318-trying-to-use-fsolve-works-over-lower-part-of-range-values-10-but-fails-in-upper-part-of-range#answer_1112058

在 MATLAB Online 中打开

 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)

5 个评论
显示 3更早的评论隐藏 3更早的评论

Walter Roberson 2022-11-26

在 MATLAB Online 中打开

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.

Walter Roberson 2022-11-26

在 MATLAB Online 中打开

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 2022-11-26

在 MATLAB Online 中打开

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 2022-11-25

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1863318-trying-to-use-fsolve-works-over-lower-part-of-range-values-10-but-fails-in-upper-part-of-range#answer_1112033

在 MATLAB Online 中打开

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 ...

3 个评论
显示 1更早的评论隐藏 1更早的评论

Paul 2022-11-27

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 2022-11-27

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.

请先登录，再进行评论。

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

更多回答（2 个）

5 个评论
显示 3更早的评论隐藏 3更早的评论

3 个评论
显示 1更早的评论隐藏 1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

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

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

更多回答（2 个）

5 个评论 显示 3更早的评论隐藏 3更早的评论

3 个评论 显示 1更早的评论隐藏 1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论

3 个评论
显示 1更早的评论隐藏 1更早的评论