Re: Q: NSolve with conditions
- To: mathgroup at smc.vnet.net
- Subject: [mg20092] Re: Q: NSolve with conditions
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Thu, 30 Sep 1999 02:43:07 -0400
- Organization: Wolfram Research, Inc.
- References: <7ssfgd$bkj@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Alphons Fakler wrote:
>
> Hello,
>
> I'm trying to solve a system of linear equations by NSolve. How can I
> restrict the solutions, to a subset, that machtes given conditions.
>
> e.g.:
> GlgSys = {EDTAtot == 0.002,
> EDTA1 == 10^10.26*EDTA*H,
> EDTA2 == 10^6.61*EDTA1*H,
> EDTA3 == 10^2.76*EDTA2*H,
> EDTA + EDTA1 + EDTA2 + EDTA3 == EDTAtot,
> H + 2*EDTAtot == 10^(-14)/H + 2*EDTA2 + 3*EDTA1 + EDTA3 + 4*EDTA};
>
> Following conditions should be kept:
> 0<= EDTA <= 0.002,
> 0<= EDTA1 <= 0.002,
> 0<= EDTA2 <= 0.002,
> 0<= EDTA3 <= 0.002,
> 0<= H <= 10^(-4)
>
> Thanks for your reply
>
> --
> Alphons Fakler
> Zentrum fur Chemische Sensoren
> ETH Technopark
> Technoparkstr. 1
> CH - 8005 Zurich
> Tel. (01) 445-1492
> Fax: (01) 445-1233
> EMail: fakler at chemsens.pharma.ethz.ch
Your first problem will be to get any solutions. This problem is
apparently quite sensitive to the large variation in scales and runs
into trouble at machine precision. To alleviate this I had success by
converting to high precision, as below. Note that I am working in
version 4 of Mathematica, offhand I do not know how earlier versions
will behave.
GlgSys = {EDTAtot - 0.002, EDTA1 - 10^10.26*EDTA*H,
EDTA2 - 10^6.61*EDTA1*H, EDTA3 - 10^2.76*EDTA2*H,
EDTA + EDTA1 + EDTA2 + EDTA3 - EDTAtot,
H + 2*EDTAtot - 10^(-14)/H + 2*EDTA2 + 3*EDTA1 + EDTA3 + 4*EDTA};
exactGlgSys = Rationalize[GlgSys,0];
vars = {EDTA,EDTA1,EDTA2,EDTA3,H,EDTAtot};
Now you can use high precision:
In[7]:= Timing[solns = NSolve[exactGlgSys, vars,
WorkingPrecision->100];]
Out[7]= {1. Second, Null}
In[8]:= Length[solns]
Out[8]= 5
Now you can get the desired solution using Select. I convert to machine
precision just to make it more concise for readability.
func[sol_] := With[{vv=vars/.sol},
0<=vv[[1]]<=.002 && 0<=vv[[2]]<=.002 &&
0<=vv[[3]]<=.002 && 0<=vv[[4]]<=.002 &&
0<=vv[[5]]<= 10^(-4)]
In[12]:= sol = N[Select[solns, func]]
-20
Out[12]= {{EDTA3 -> 4.89993 10 , EDTAtot -> 0.002, EDTA -> 0.00197005,
-10 -13
> EDTA1 -> 0.0000299489, EDTA2 -> 1.01926 10 , H -> 8.35418 10 }}
There is an undocumented experimental technique that allows NSolve to
weed out partial or full candidate solutions as it produces them. This
was developed not so much for convenience but rather for efficiency; In
some instances, say when one seeks only real solutions to a large
polynomial system, the time savings can be tremendous. This is not
really an issue here, where it takes a second to get all solutions, and
essentially no time to discard the unwanted ones. That said, here goes.
Were it not hidden, the usage message would be something like this:
SelectCriterion::usage =
"SelectCriterion is an option for NSolve that determines whether to
keep or discard a root in any variable. The value provided should be a
pure
function. It will be applied to the pair {var, candidate value} to
determine
whether to keep or discard that value. Hence one can use different
criteria
for different variables, provided these criteria are independent of one
another."
To use it in this problem:
In[16]:= soltoo = N[NSolve[N[exactGlgSys,100], vars,
WorkingPrecision->100,
Internal`SelectCriterion->(Switch[#[[1]],
EDTA || EDTA1 || EDTA2 || EDTA3, 0<=#[[2]]<=0.002,
H, 0<=#[[2]]<=10^(-4),
_, True]&)]]
-20
Out[16]= {{EDTA3 -> 4.89993 10 , EDTAtot -> 0.002, EDTA -> 0.00197005,
-10 -13
> EDTA1 -> 0.0000299489, EDTA2 -> 1.01926 10 , H -> 8.35418 10 }}
In[17]:= soltoo===sol
Out[17]= True
Obviously such undocumented features are subject to change or removal;
this should be kept in mind if you intend to use them.
Daniel Lichtblau
Wolfram Research