4
$\begingroup$

I want to ask how to do the numerical double integral like this (Actually, I want to carry out the parametric plot of Trial1[t] and Trial2[t] as t varies. But, since Trial2[tt] cannot be calculated, I just consider Trial2[2.5] first. )

Trial1[t_?NumericQ] := Exp[NIntegrate[Log[Log[x]], {x, 2, t}]];

Trial2[tt_?NumericQ] := Integrate[Exp[Exp[tt]] Trial1[tt], {tt, 2, rr}];

Trial2[2.5] // N

Another follow-up question is that if my code involves NDsolve like the following. (Sorry for very long code.)

tpa = 6;
a1 = (tpa + 1.08) Pi/(tpa + 0.1);
b1 = (tpa + 1.8) Pi/(tpa + 1.3);
a2 = tpa Pi/(tpa + 1.19);
b2 = (tpa + 1.37) Pi/(tpa + 2);
m1 = 0.45;
m2 = 0.55;
lambda1 = 2.6*(10^(-4));
lambda2 = lambda1*(m1/m2)^2;
Mpl = 1;

(*Define the saxion. *)
chi10 =  (b2 Log[lambda1] - b1 Log[lambda2])/(a2 b1 - a1 b2) ;
chi20 = - (a2 Log[lambda1] - a1 Log[lambda2])/(a2 b1 - a1 b2) ;


(*Define the decay constants. *)
f1 = 1/a1 ;
f2 = 1/b1 ;
g1 = 1/a2 ;
g2 = 1/b2 ;

(*Define the Amplitudes. *)
Lambda1 = ((m1^4 lambda1^2)/(2 k chi10 chi20))^(1/4) ;
Lambda2 = ((m2^4 lambda2^2)/(2 k chi10 chi20))^(1/4) ;


(*Effective decay constants after KNP mechanism *)
fpsi = (f1 f2)/(Sqrt[(f1)^2 + (f2)^2]) ;
fpsiprime = (g1 g2 (Sqrt[(f1)^2 + (f2)^2]))/(f1 g1 + f2 g2) ;
feff = (g1 g2 (Sqrt[(f1)^2 + (f2)^2]))/(Abs[f1 g2 - f2 g1]) ;

(*Define the parameters for the calculations.*)
(*Potential Derivatives*)
V[phi1_, phi2_] := (Lambda1)^(4) ( 1 - Cos[ phi1/ fpsi] ) + (Lambda2)^(4) ( 1 - Cos[ phi1/fpsiprime + phi2/feff ] ) ;

(*1st Derivatives*)
V1[phi1_, phi2_] = D[V[phi1, phi2], {phi1, 1}];
V2[phi1_, phi2_] = D[V[phi1, phi2], {phi2, 1}];
(*2nd Derivatives*)
V11[phi1_, phi2_] = D[V[phi1, phi2], {phi1, 2}];
V12[phi1_, phi2_] = D[V[phi1, phi2], phi1, phi2];
V22[phi1_, phi2_] = D[V[phi1, phi2], {phi2, 2}];

(*Define the rotation angle.*)
Theta[phi1_, phi2_] := ArcTan[V1[phi1, phi2], V2[phi1, phi2]];

(*Define the rotations of the 1st derivatives of potential.*)
VSigma[phi1_, phi2_] := V1[phi1, phi2]*Cos[Theta[phi1, phi2]] + V2[phi1, phi2]*Sin[Theta[phi1, phi2]];

(*Define the slow-roll parameters.*)
(*1st order slow-roll parameters.*)
Epsilon11[phi1_, phi2_] := (1/2) (Mpl)^2 (V1[phi1, phi2]/V[phi1, phi2])^2 ; 
Epsilon12[phi1_, phi2_] := (1/2) (Mpl)^2 (V1[phi1, phi2]/V[phi1, phi2]) (V2[phi1, phi2]/V[phi1, phi2]) ;
Epsilon22[phi1_, phi2_] := (1/2) (Mpl)^2 (V2[phi1, phi2]/V[phi1, phi2])^2 ;

(*2nd order slow-roll parameters.*)
Eta11[phi1_, phi2_] := (Mpl)^2 V11[phi1, phi2]/V[phi1, phi2] ;
Eta12[phi1_, phi2_] := (Mpl)^2 V12[phi1, phi2]/V[phi1, phi2] ;
Eta22[phi1_, phi2_] := (Mpl)^2 V22[phi1, phi2]/V[phi1, phi2] ;

(*High order derivative terms of potential.*)
(*2nd order*)
VSigmaSigma [phi1_, phi2_] := (Cos[Theta[phi1, phi2]]^2) V11[phi1, phi2] + 2 (Sin[Theta[phi1, phi2]]) (Cos[Theta[phi1, phi2]]) V12[phi1, phi2] + (Sin[Theta[phi1, phi2]]^2) V22[phi1, phi2] ;
VSigmaS[phi1_, phi2_] := - Sin[Theta[phi1, phi2]]Cos[Theta[phi1,phi2]]V11[phi1, phi2] + (Cos[Theta[phi1, phi2]]^2 -Sin[Theta[phi1, phi2]]^2) V12[phi1, phi2] + Sin[Theta[phi1, phi2]]Cos[Theta[phi1, phi2]] V22[phi1, phi2] ;
VSS[phi1_, phi2_] := Sin[Theta[phi1, phi2]]^2 V11[phi1, phi2] - 2Sin[Theta[phi1, phi2]] Cos[Theta[phi1, phi2]] V12[phi1, phi2] + Cos[Theta[phi1, phi2]]^2 V22[phi1, phi2] ;

(*Slow-roll parameters after the rotations to sigma and s.*)
(*1st order slow-roll parameters.*)
Epsilon[phi1_, phi2_] := Epsilon11[phi1, phi2] + Epsilon22[phi1, phi2];

(*2nd order slow-roll parameters.*)
EtaSigmaSigma[phi1_, phi2_] := Cos[Theta[phi1, phi2]]^2 Eta11[phi1, phi2] + 2 Sin[Theta[phi1, phi2]] Cos[Theta[phi1, phi2]] Eta12[phi1, phi2] + Sin[Theta[phi1, phi2]]^2 Eta22[phi1, phi2] ;
EtaSigmaS[phi1_, phi2_] := - Sin[Theta[phi1, phi2]] Cos[Theta[phi1, phi2]] Eta11[phi1, phi2] + (Cos[Theta[phi1, phi2]]^2 - Sin[Theta[phi1, phi2]]^2) Eta12[phi1, phi2] + Sin[Theta[phi1, phi2]] Cos[Theta[phi1, phi2]] Eta22[phi1, phi2] ;
EtaSS[phi1_, phi2_] := Sin[Theta[phi1, phi2]]^2 Eta11[phi1, phi2] - 2 Sin[Theta[phi1, phi2]] Cos[Theta[phi1, phi2]] Eta12[phi1, phi2] + Cos[Theta[phi1, phi2]]^2 Eta22[phi1, phi2] ;

(*Numerical calculations of phi1 and phi2*)
ww = - 0.95;  (* As required, w < -1/3 *)
numsol = NDSolve[{phi1''[t] + 3(2/(3(1+ww)t))phi1'[t] +((Lambda1)^4/fpsi) Sin[(phi1[t]/fpsi)] == 0, phi2''[t] + 3(2/(3(1+ww)t))phi2'[t] + ((Lambda2)^4/fpsiprime) Sin[(phi1[t]/fpsiprime) +(phi2[t]/feff)] == 0, phi1[1] == 4.374641, phi2[1] == 15.2, phi1'[1] == 0, phi2'[1] == 0}, {phi1, phi2}, {t, 0.01, 100000}] ; (* phi'’+3 H phi’ + V’ == 0 *) 

Alpha[t_] := 2 EtaSigmaS[phi1, phi2] /. {phi1 -> Evaluate[phi1[t] /. numsol], phi2 -> Evaluate[phi2[t] /. numsol]} ;
BBeta[t_] := -2 Epsilon[phi1, phi2] + EtaSigmaSigma[phi1, phi2] - EtaSS[phi1, phi2] /. {phi1 -> Evaluate[phi1[t] /. numsol], phi2 -> Evaluate[phi2[t] /. numsol]} ;
TSS[t_ /; t > 0.999] = Exp[Integrate[BBeta[x] (2/(3 (1 + ww) x)), {x, 1, t},Assumptions -> {t > 0.999}]];
TRS[t_?NumericQ] := Integrate[Alpha[x] (2/(3 (1 + ww) x)) TSS[x], {x, 1, t}];
ParametricPlot[{TSS[t], TRS[t]}, {t, 1, 1.1}, AspectRatio -> 1]

In the last line, when I want to do the parametric plot of TSS[t] and TRS[t] faster as t varies, how should I change the code in TSS[t] and TRS[t]? Or I need to change the code in other places ?

Thank you very much for your help.

$\endgroup$
  • $\begingroup$ Your code for Trial2 is not correct. Did you mean Trial2[t_?NumericQ] := Integrate[Exp[Exp[x]] Trial1[x], {x, 2, t}]? $\endgroup$ – KraZug Jun 15 '17 at 11:34
  • $\begingroup$ Yes, in particular, I want to set up the function as Trial[t_?NumericQ]:=Integrate[Exp[Exp[x]] Trial1[x], {x, 2, t}]. But, I am afraid that the computer will misunderstand that the variable t here is the same as that in Trial1 function. So, I change it into variable tt here instead of variable t. (I am a beginner of mathematica. ) $\endgroup$ – Ellgan Jun 15 '17 at 15:25
  • $\begingroup$ In your updated code, you need to learn the difference between := (SetDelayed) and = (Set). In short, SetDelayed will only evaluate the RHS when it is called later on, while Set will evaluate it the RHS once and keep that value no matter if the RHS is later changed. $\endgroup$ – KraZug Jun 15 '17 at 19:21
6
$\begingroup$

As Jack LaVigne pointed out, Mathematica can do the first integral analytically for $t>1$. Which means you can write (with Set instead of SetDelayed)

Trial1[t_/;t>1] = Exp[Integrate[Log[Log[x]], {x, 2, t}, Assumptions -> {t > 1}]]; 
Trial2[t_?NumericQ] := NIntegrate[Exp[Exp[x]] Trial1[x], {x, 2, t}];

Which will calculate Trial1[t] once and return an expression which is valid for all $t>1$. This speeds up the calculation of Trial2[t] hugely, and it now takes 0.01 second per value of t on my machine.

For $t>2$, Trial2 gets doubly exponentially large, but you can plot it with ParametricPlot easily enough:

ParametricPlot[{Trial1[t], Log@Trial2[t]}, {t, 2.1, 3}, AspectRatio -> 1] // Timing

(* 3.7 seconds *)

enter image description here

$\endgroup$
  • 1
    $\begingroup$ Good answer. You can make it run faster by using Integrate rather than NIntegrate for Trial1 and changing t_?NumericQ to t_. It gives a closed from answer provided that t>1. Suggest using the definition Trial1[t_ /; t > 1] := Exp[-2 Log[Log[2]] + t Log[Log[t]] + LogIntegral[2] - LogIntegral[t]] $\endgroup$ – Jack LaVigne Jun 15 '17 at 14:10
  • $\begingroup$ @JackLaVigne, I didn't realise an analytic result existed. $\endgroup$ – KraZug Jun 15 '17 at 15:48
  • $\begingroup$ @JackLaVigne, thanks, I've incorporated that fact into my answer. $\endgroup$ – KraZug Jun 15 '17 at 15:57

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.