Message Boards Message Boards

Producing a smooth curve from numerical minimization.

GROUPS:

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)/(2lA) - (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.

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 to this post.

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]
Attachments:
POSTED BY: Deb Chicago
Answer
9 days ago

Use

WorkingPrecision -> 30

added: To the FindMinimum[]

Change the constants to have consistent precision:

phiAlist = Range[0.01`30, 0.99`30, 0.01`30];
totalpoints = Length[phiAlist];
chi = 0.0004`30;
epsilonA = 0;
epsilonB = 0;
epsilonAB = 0;

and use ListLinePlot instead of ListPlot and they are smooth. Although I do not understand how you get the piecewise function from the original function.

Regards

POSTED BY: Neil Singer
Answer
8 days ago

Group Abstract Group Abstract