0
$\begingroup$

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;
cc1R[R_]:=cc1((1.0+a(cp*R/cosi^2)^2+b(cp*R/cosi^2)^3)/(1.0+c(cp*R/cosi^2)^2+d(cp*R/cosi^2)^4));

S0=0.210; a0=3.51; a1=-2.47; a2=15.2; b3=-1.75;
aS[S_]:=a0+a1(S0/S)+a2(S0/S)^2;
bS[S_]:=b3(S0/S)^3;

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*);

(*Integrand-I*)
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];

(*Integrand-II*)
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];

(*Integrand-III*)
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.

Improve this question
$\endgroup$
7
  • $\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$
    – Chowdhury Aminul Islam
    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$
    – Chowdhury Aminul Islam
    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

0

Reset to default

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.