The parameters. Here I mention all the parameters I use in my code including some of the functions that I define for convenience

Parameter set-I

my=0.0055; my3=0.1407; co=0.6023; cc1=3.67/co^2; cc2=12.36/co^5;

cosi=0.3; a=0.0108805; b=-1.0133*10^-4; c=0.02228; d=1.84558*10^-4;

S0=0.210; a0=3.51; a1=-2.47; a2=15.2; b3=-1.75;

bhpi=0.140; cp=0.302862; cp1=2.0/3.0 cp; cp2=1.0/3.0 cp;

Predefined functions set-I

My1[y1_, y2_, y3_, R_]:=my-2.0*cc1R[R]*y1+2.0*cc2*y2*y3; 
My2[y1_, y2_, y3_, R_]:=my-2.0*cc1R[R]*y2+2.0*cc2*y1*y3;
My3[y1_, y2_, y3_, R_]:=my3-2.0*cc1R[R]*y3+2.0*cc2*y1*y2;
DPy1[bb_, y1_, y2_, y3_, R_]:=(bb^2+My1[y1, y2, y3, R]^2)^(1/2)
DPy2[bb_, y1_, y2_, y3_, R_]:=(bb^2 + My2[y1, y2, y3, R]^2)^(1/2)
DPy3[bb_, y1_, y2_, y3_, R_]:=(bb^2 + My3[y1, y2, y3, R]^2)^(1/2)

DPy1Rcc1R[bbz_, y1_, y2_, y3_, R_, l_]:=(bbz^2+My1[y1, y2, y3, R]^2+2.0*l cp1 R)^(1/2);
DPy2Rcc1R[bbz_, y1_, y2_, y3_, R_, l_]:=(bbz^2+My2[y1, y2, y3, R]^2+2.0*l cp2 R)^(1/2);
DPy3Rcc1R[bbz_, y1_, y2_, y3_, R_, l_]:=(bbz^2+My3[y1, y2, y3, R]^2+2.0*l cp2 R)^(1/2);

xy1cc1R[y1_, y2_, y3_, R_]:=My1[y1, y2, y3, R]^2/(2.*cp1*R);
xy2cc1R[y1_, y2_, y3_, R_]:=My2[y1, y2, y3, R]^2/(2.*cp2*R);
xy3cc1R[y1_, y2_, y3_, R_]:=My3[y1, y2, y3, R]^2/(2.*cp2*R);

The definition of the parameters and functions ends here. Bellow I give the expressions for the final equation and its different parts.

(*Part I of the final function*)
plp[y4_, S_]:=S^4(-(aS[S]/2.0)y4^2+bS[S]Log[1-6.0 y4^2+8.0 y4^3-3.0 y4^4])(*cannot directly put y4=1 because of the log term*);

int1[y1_?NumericQ, y2_?NumericQ, y3_?NumericQ, y4_?NumericQ, S_?NumericQ, R_?NumericQ, tp_?NumericQ, l_?IntegerQ]:=NIntegrate[(1.0/(1.0-tp)) ((1.0+3.0 y4(1.0+(1.+(1.-tp)(-(DPy1Rcc1R[bbz, y1, y2, y3, R, l]/S)))^(1./(1.-tp)))(1.+(1.-tp)(-(DPy1Rcc1R[bbz, y1, y2, y3, R, l]/S)))^(1./(1.-tp))+(1.+(1.-tp)(-3.0 DPy1Rcc1R[bbz, y1, y2, y3, R, l]/S))^(1./(1.- tp)))^(1.0-tp)-1.0), {bbz, -\[Infinity], \[Infinity]},Method -> {Automatic,"SymbolicProcessing" -> False},AccuracyGoal -> 10];

int2[y1_?NumericQ, y2_?NumericQ, y3_?NumericQ, y4_?NumericQ, S_?NumericQ, R_?NumericQ, tp_?NumericQ,l_?IntegerQ]:=NIntegrate[(1.0/(1.0-tp))((1.0+3.0 y4 (1.0+(1.+(1.-tp)(-(DPy2Rcc1R[bbz, y1, y2, y3, R, l]/S)))^(1./(1.-tp)))(1.+(1.-tp)(-(DPy2Rcc1R[bbz, y1, y2, y3, R, l]/S)))^(1./(1.-tp))+(1.+(1.-tp)(-3.0 DPy2Rcc1R[bbz, y1, y2, y3, R, l]/S))^(1./(1.-tp)))^(1.0 - tp)-1.0),{bbz, -\[Infinity], \[Infinity]},Method -> {Automatic,"SymbolicProcessing" -> False},AccuracyGoal -> 10];

int3[y1_?NumericQ, y2_?NumericQ, y3_?NumericQ, y4_?NumericQ, S_?NumericQ, R_?NumericQ, tp_?NumericQ,l_?IntegerQ]:=NIntegrate[(1.0/(1.0-tp))((1.0+3.0 y4(1+(1.+(1.-tp)(-(DPy3Rcc1R[bbz, y1, y2, y3, R, l]/S)))^(1./(1.-tp)))(1.+(1.-tp)(-(DPy3Rcc1R[bbz, y1, y2, y3, R, l]/S)))^(1./(1.-tp))+(1.+(1.-tp)(-3.0 DPy3Rcc1R[bbz, y1, y2, y3, R, l]/S))^(1./(1.-tp)))^(1.0-tp)-1.0),{bbz, -\[Infinity], \[Infinity]},Method -> {Automatic,"SymbolicProcessing" -> False}, AccuracyGoal -> 10];

(*The final function of which I try to find the minimum is written bellow*)
fe[y1_?NumericQ, y2_?NumericQ, y3_?NumericQ, y4_?NumericQ, 
  S_?NumericQ, R_?NumericQ, tp_?NumericQ, nl_?IntegerQ]:=cc1R[R] (y1^2 + y2^2 + y3^2)-4.cc2*y1*y2*y3-3.0 /\[Pi]^2 NIntegrate[bb^2*DPy1[bb, y1, y2, y3, R], {bb, 0, co},     Method -> {Automatic, "SymbolicProcessing" -> False}, AccuracyGoal -> 10]-3.0 /\[Pi]^2 NIntegrate[bb^2*DPy2[bb, y1, y2, y3, R], {bb, 0, co},Method -> {Automatic,"SymbolicProcessing" -> False},AccuracyGoal -> 10]-3.0 /\[Pi]^2 NIntegrate[bb^2*DPy3[bb, y1, y2, y3, R], {bb, 0, co},Method -> {Automatic, "SymbolicProcessing" -> False}, AccuracyGoal -> 10]-3.0/(2. \[Pi]^2)(cp1*R)^2(xy1cc1R[y1, y2, y3, R]/2.(3./2. xy1cc1R[y1, y2, y3, R]-1)+ xy1cc1R[y1, y2, y3, R]/2.Log[xy1cc1R[y1, y2, y3, R]/(2.0 \[Pi])]-xy1cc1R[y1, y2, y3, R]^2/2. Log[xy1cc1R[y1, y2, y3, R]]+PolyGamma[-2, xy1cc1R[y1, y2, y3, R]])-3.0/(2. \[Pi]^2)(cp2*R)^2(xy2cc1R[y1, y2, y3, R]/2.(3./2. xy2cc1R[y1, y2, y3, R]-1)+xy2cc1R[y1, y2, y3, R]/ 2.Log[xy2cc1R[y1, y2, y3, R]/(2.0 \[Pi])] -xy2cc1R[y1, y2, y3, R]^2/2. Log[xy2cc1R[y1, y2, y3, R]]+PolyGamma[-2, xy2cc1R[y1, y2, y3, R]])-3.0/(2. \[Pi]^2)(cp2*R)^2(xy3cc1R[y1, y2, y3, R]/2. (3./2. xy3cc1R[y1, y2, y3, R]-1)+      xy3cc1R[y1, y2, y3, R]/2.Log[xy3cc1R[y1, y2, y3, R]/(2.0 \[Pi])]-xy3cc1R[y1, y2, y3,R]^2/2. Log[xy3cc1R[y1, y2, y3, R]] +PolyGamma[-2, xy3cc1R[y1, y2, y3, R]])-1.0/(2. \[Pi]^2) S*cp1*R*Sum[ (2 - KroneckerDelta[l, 0]) int1[y1, y2, y3, y4, S, R, tp,l], {l, 0, nl}]-   1.0/(2. \[Pi]^2) S*cp2*R*Sum[ (2 - KroneckerDelta[l, 0]) int2[y1, y2, y3, y4, S, R, tp,       l], {l, 0, nl}]-1.0/(2. \[Pi]^2) S*cp2*R*Sum[ (2 - KroneckerDelta[l, 0]) int3[y1, y2, y3, y4, S, R, tp,l], {l, 0, nl}] + plp[y4, S]
(*The definition of function ends here.*)

Now we try to find out the minimum for the four different variables, y1, y2, y3 and y4.

y1FM = -0.001; y2FM = -0.001; y3FM = -0.001; y4FM = 0.9;(*These are the initial choices*)
tab = ParallelTable[{S, {y1, y2, y3, y4} /.FindMinimum[fe[y1, y2, y3, y4, S, 40 bhpi^2/cp, 1.1,20], {{y1, y1FM}, {y2, y2FM}, {y3, y3FM}, {y4, y4FM}}, Method -> "InteriorPoint"][[2]]}, {S, 0.0001, 0.4, 0.001}]

When I generate the table by finding the minimum of the function 'fe' through the use of FindMinium, all the solutions (y1, y2, y3, y4) possess a lot of unphysical oscillations! enter image description here

y2, y3 and y4 also possess similar oscillatory points. How can I get rid of those? Please let me know. Any help or clue is much appreciated. Thanks in advance.

  • $\begingroup$ Running your code, the first discontinuity appears between: {0.0081, {-0.0268738, -0.0174615, -0.0193506, -0.000192137} and {0.0091, {-1.51743, 1.18273, -1.00136, 6.60748*10^-6}. Looks to me as if nothing has gone wrong, it's the nature of your function that the minima jump around. $\endgroup$
    – bill s
    Mar 18 at 21:17
  • $\begingroup$ Thank you very much for your comment. Actually, the function I work with is a physically motivated function. It should in principle produce a result like the solid line shown in the figure I shared, except for those oscillatory points. Also changing the method changes those oscillatory points' numbers and positions. Thus, I wonder whether I am doing something wrong with the numerical treatment. $\endgroup$ Mar 19 at 12:28
  • $\begingroup$ I would look carefully at the terms that raise variables to a fractional power. Since powers need not be unique (simple example, sqrt has 2 answers) it seems plausible that maybe the FindMinimum may be switching between the two values. $\endgroup$
    – bill s
    Mar 19 at 13:41
  • $\begingroup$ I checked further with the function. For some of the values of the variables it becomes complex. I don't know whether that is causing the problem. The function should in principle be always real. I tried to constrain the evaluation by restricting the function to be always real. Though it runs and gives expected result for a smaller number of points, but as I increase the sampling points the programme keeps on running. $\endgroup$ Mar 21 at 11:22
  • $\begingroup$ If the answers are becoming complex, maybe it can be fixed by wrapping those terms in Abs[ ] thereby forcing it to be real? Alternatively, if the imaginary parts are small, you could remove them with Chop[ ]. $\endgroup$
    – bill s
    Mar 21 at 13:27


Your Answer

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

Browse other questions tagged or ask your own question.