Numerically solving a system of nonlinear equations - chemical equilibrium concentrations
5 views (last 30 days)
Show older comments
Hi-
Can anyone help me fix this?
KMD = 1.3e-3; % 2 monomer -> dimer
KMT = 4.1e-12; % 4 monomer -> tetramer
KTO = 2.4e-17; % 4 monomer + tetramer -> octamer
ctot = 0.005; % total concentration (monomer)
func=@(C) ([ (C(1)^2 / C(2))/KMD - 1;...
(C(1)^4 / C(3))/KMT - 1; ...
(C(3) * C(1)^4 /C(4))/KTO - 1; ...
(C(1)+C(2)*2+C(3)*4+C(4)*8)/ctot-1]) ;
x0 = [1/4 1/8 1/16 1/32]' *ctot;
[xpar,fehler,resid,flag] = ...
fminunc(@(x) sum(func(x)),x0);
xpar*1e3
fehler
I get this error:
Warning: Gradient must be provided for trust-region method;
using line-search method instead.
> In fminunc at 265
Line search cannot find an acceptable point along the current
search direction.
Thanks for any help!
0 Comments
Answers (1)
Walter Roberson
on 29 Aug 2011
I notice that you are using fminunc() for a system that involves concentrations. I would have thought that you would not like any concentration to go below 0 or above 1, and thus would be more likely to want to use fmincon ?
In the system as-is, unconstrained, the system has a minimum of -infinity when C2=-infinity, or when C4=-infinity, or when C3=infinity and C4 is sufficiently negative, or when C3=-infinity and C4 is positive or insufficiently negative... amongst other conditions.
If the concentrations are constrained to 0 at minimum, then the system is minimized when C1=0 and C2, C3, and C4 approach 0 at the limit (but are non-zero). Another way of putting that is that although the system is sigular if C2 or C3 or C4 is 0, it is a "removable singularity" if they are all 0, with the equation becoming -4 + 0/0
2 Comments
Walter Roberson
on 30 Aug 2011
Add the Option to check the function values. With you using 0 as your initial concentrations, you are going to get nan out of your func, and it isn't going to be able to figure out where to go.
Also, your func only outputs a column of 3 values: you lost the output about the total concentrations. That will cause problems.
I would suggest that you do not want 0 as lb because of the division by 0 those values get. Use a small non-negative value instead, and be sure that your x0 falls within those constraints.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!