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

1 次查看(过去 30 天)
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)

采纳的回答

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

更多回答(2 个)

Walter Roberson
Walter Roberson 2022-11-25
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 个评论
Walter Roberson
Walter Roberson 2022-11-26
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
Walter Roberson 2022-11-26
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.)

请先登录,再进行评论。


Paul
Paul 2022-11-25
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 个评论
Paul
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
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.

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Solver Outputs and Iterative Display 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by