0
$\begingroup$

I have a function f(pA, pB, pAB, x) where pA,pB and pAB depend on x.

f[pA, pB, pAB, x] = -(x pA/(2 lA)) Log[(8 x)/(lA E)] - ((1 - x) pB/(2 lB)) Log[8 (1 - x)/(lB E)] - (x pAB/ lA) Log[8(1 - x)/(lB E)] + (x/(2 lA)) (pA Log[pA] + 2 pAB Log[AB] + 2 (1 - pA - pAB) Log[1 - pA - pAB]) + ((1 - x)/(2 lB)) (pB Log[pB] + 2(1 - pB - ((lB x)/(lA (1 - x))) pAB) Log[(1 - pB - ((lB x)/(lA (1 - x))) pAB)]) - (epsilonA x pA)/(2 lA) - (epsilonB (1 - x) pB)/(2 lB) - (epsilonAB x pAB)/lA

pA,pB and pAB are found by minimizing f at each x with respect to pA,pB and pAB.

However, there are constraints such as pA<1, pB<1, pAB<1, pA+pAB<1 and pB+pAB*x/(1-x)<=1.

lA, lB, epsilonA, epsilonB, epsilonAB are constants.

My goal is to plot f as a function of x in the range (0,1).

I used FindMinimum to minimize f with respect to pA,pB and pAB at each x (in a loop) and then plug the solution (pA, pB and pAB) into f to get its value. However, this leads to curves with kinks. This is the last thing I would like to see. I used PrecisionGoal, SetPrecision, but to no avail.The file is attached here. I would also like to bring to your notice that when I perform constrained minimization using FindMinimum does not obey the constraints while exploring the search space, (hence resulting in errors such as the function is not real for ..). This caused me to resort to using a piecewise function, where I use the limiting values of the function at the boundary such as for pA+pAB<1, etc.. The code is attached below if you would like to copy it and test it. Note that in the code, phiAloc is the same as x that I defined in f [pA,pB,pAB,x].

lA = 500;
lB = 500;

NA = 5000;
NB = 5000;
phiAlist = Range[0.01, 0.99, 0.01];
totalpoints = Length[phiAlist];
chi = 0.0004;
epsilonA = 0.0;
epsilonB = 0.0;
epsilonAB = 0.0;
funcfint[pA_, pB_, pAB_, 
phiAloc_] := (Piecewise[{{-(phiAloc*pA/(2*lA))*
   Log[(8*phiAloc)/(lA*E)] - ((1 - phiAloc)*pB/(2*lB))*
   Log[8*(1 - phiAloc)/(lB*E)] - (phiAloc*pAB/lA)*
   Log[8*(1 - phiAloc)/(lB*E)] + (phiAloc/(2*lA))*(pA*Log[pA] + 
     2 pAB*Log[pAB] + 
     2*(1 - pA - pAB)*
      Log[1 - pA - pAB]) + ((1 - phiAloc)/(2*lB))*(pB*Log[pB] + 
     2*(1 - pB - ((lB*phiAloc)/(lA*(1 - phiAloc)))*pAB)*
      Log[(1 - 
         pB - ((lB*phiAloc)/(lA*(1 - phiAloc)))*
          pAB)]) - (epsilonA*phiAloc*pA)/(2*
     lA) - (epsilonB*(1 - phiAloc)*pB)/(2*lB) - (epsilonAB*
     phiAloc*pAB)/lA, 
 pA + pAB < 1 && 
  pB + pAB*(lB * phiAloc)/(lA*(1 - phiAloc)) < 
   1}, { -(phiAloc*pA/(2*lA))*
   Log[(8*phiAloc)/(lA*E)] - ((1 - phiAloc)*pB/(2*lB))*
   Log[8*(1 - phiAloc)/(lB*E)] - (phiAloc*pAB/lA)*
   Log[8*(1 - phiAloc)/(lB*E)] + (phiAloc/(2*lA))*(pA*Log[pA] + 
     2 pAB*Log[pAB] + 0) + ((1 - phiAloc)/(2*lB))*(pB*Log[pB] + 
     2*(1 - pB - ((lB*phiAloc)/(lA*(1 - phiAloc)))*pAB)*
      Log[(1 - 
         pB - ((lB*phiAloc)/(lA*(1 - phiAloc)))*
          pAB)]) - (epsilonA*phiAloc*pA)/(2*
     lA) - (epsilonB*(1 - phiAloc)*pB)/(2*lB) - (epsilonAB*
     phiAloc*pAB)/lA, 
 pA + pAB >= 
   1 && (pB + pAB*(lB * phiAloc)/(lA*(1 - phiAloc)) < 
    1)}, { -(phiAloc*pA/(2*lA))*
   Log[(8*phiAloc)/(lA*E)] - ((1 - phiAloc)*pB/(2*lB))*
   Log[8*(1 - phiAloc)/(lB*E)] - (phiAloc*pAB/lA)*
   Log[8*(1 - phiAloc)/(lB*E)] + (phiAloc/(2*lA))*(pA*Log[pA] + 
     2 pAB*Log[pAB] + 
     2*(1 - pA - pAB)*
      Log[1 - pA - pAB]) + ((1 - phiAloc)/(2*lB))*(pB*Log[pB] + 
     0) - (epsilonA*phiAloc*pA)/(2*lA) - (epsilonB*(1 - phiAloc)*
     pB)/(2*lB) - (epsilonAB*phiAloc*pAB)/lA, 
 pA + pAB < 1 && 
  pB + pAB*(lB * phiAloc)/(lA*(1 - phiAloc)) >= 
   1}, { -(phiAloc*pA/(2*lA))*
   Log[(8*phiAloc)/(lA*E)] - ((1 - phiAloc)*pB/(2*lB))*
   Log[8*(1 - phiAloc)/(lB*E)] - (phiAloc*pAB/lA)*
   Log[
    8*(1 - phiAloc)/(lB*E)] + (phiAloc/(2*lA))*(pA*Log[pA] + 
     2 pAB*Log[pAB] + 0) + ((1 - phiAloc)/(2*lB))*(pB*Log[pB] + 
     0) - (epsilonA*phiAloc*pA)/(2*lA) - (epsilonB*(1 - phiAloc)*
     pB)/(2*lB) - (epsilonAB*phiAloc*pAB)/lA, 
 pA + pAB >= 1 && 
  pB + pAB*(lB * phiAloc)/(lA*(1 - phiAloc)) >= 1}}])
finterlist = {};
composition = {};
pABlist = {};
pAlist = {};
pBlist = {};
sumlist = {};
Do[ phiAloc = phiAlist[[i]]; 
  g = FindMinimum[{funcfint[u, v, w, phiAloc], 
 0 < u <= 1 && 0 < v <= 1 && 0 < w <= 1 && u + w <= 1 && 
  v + (lB*phiAloc/(lA*(1 - phiAloc)))*w <= 1}, {u, v, w}] ;
  prob = {u, v, w} /. g[[2]];
  pA = prob[[1]];
  pB = prob[[2]];
  pAB = prob[[3]];
  sum = pB + (phiAloc/(1 - phiAloc))*
 pAB;(*just checking whether the constraint is satisfied*)

  finter = funcfint[pA, pB, pAB, phiAloc];
  (* free=funcfenergy[pA,pB,pAB,phiAloc];*)

  AppendTo[composition, phiAloc]; AppendTo[pAlist, pA]; 
  AppendTo[pBlist, pB]; AppendTo[pABlist, pAB]; 
  AppendTo[finterlist, finter];
  AppendTo[sumlist, sum],
  {i, 1, totalpoints, 1}];
finterplot = Transpose[{composition, finterlist}];
finterfunction = Interpolation[finterplot, Method -> "Spline"];
pAplot = Transpose[{composition, pAlist}];
pBplot = Transpose[{composition, pBlist}];
pABplot = Transpose[{composition, pABlist}];
sumplot = Transpose[{composition, sumlist}];
Plot[finterfunction[x], {x, 0.01, 0.99}]
Plot[finterfunction'[x], {x, 0.01, 0.99}]
ListPlot[pAplot]
ListPlot[pBplot]
ListPlot[pABplot]
ListPlot[sumplot]
$\endgroup$
  • 1
    $\begingroup$ It is not likely that you will receive many responses with such a long block of code. Can you illustrate your issue with a shorter example? $\endgroup$ – bbgodfrey Jul 14 '17 at 22:13
  • $\begingroup$ Cross-posted: community.wolfram.com/groups/-/m/t/1146934 $\endgroup$ – Michael E2 Jul 15 '17 at 2:39
1
$\begingroup$

The problem was solved on http://community.wolfram.com/groups/-/m/t/1146934 by setting WorkingPrecision -> 30 and setting the precision of the parameters to 30, too.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.