Numerically solving a system of nonlinear equations - chemical equilibrium concentrations

5 views (last 30 days)
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!

Answers (1)

Walter Roberson
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
Martin Allan
Martin Allan on 30 Aug 2011
Thank you Walter!
Of course, concentrations should be positive.
Here's an updated version using fmincon:
clear all
close all
KMD = 1.3e-3; % 2 monomer -> dimer
KMT = 4.1e-12; % 4 monomer -> tetramer
KTO = 2.4e-17; % 4 monomer + tetramer -> octamer
ctot = 0.0001; % 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]) ;
% Total concentration == ctot:
A = [ 1 2 4 8;...
-1 -2 -4 -8];
b=[ctot ctot]';
Aeq=[]; beq=[];
% no negative concentrations:
lb = [0 0 0 0]';
% This should be redundant:
ub = [1 1/2 1/4 1/8]' * ctot;
% ub = [Inf Inf Inf Inf]';
nonlcon = [];
% x0 = [1/4 1/8 1/16 1/32]' *ctot;
x0 = [ctot 0 0 0]';
options = optimset('MaxFunEvals',1e5,...
'MaxIter',500);
x = fmincon(@(x) sum(func(x).^2),x0,A,b,Aeq,beq,lb,ub,nonlcon,options)
ctot2 = [1 2 4 8] * x
fval = [ (x(1)^2 / x(2))/KMD;...
(x(1)^4 / x(3))/KMT; ...
(x(3) * x(1)^4 /x(4))/KTO]
Unfortunately, this does not work as I would like it to.
Warning: Large-scale (trust region) method does not currently solve this type of problem,
using medium-scale (line search) instead.
is probably OK;
> In fmincon at 317
Maximum number of iterations exceeded;
increase OPTIONS.MaxIter.
is more worrying. I tried increasing MaxIter up to 50,000, but the result is always:
x =
1.0e-04 *
1.0000
0
0
0
which does not fit my criteria.
Any thoughts?
Walter Roberson
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.

Sign in to comment.

Categories

Find more on Chemistry in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!