0
$\begingroup$

I am doing this integration and when I use $C_\alpha(T)$ as just constant numbers(Subscript[C, \[Alpha]][T_]= Piecewise[{{1.2, T < 100}, {8, T > 1000}, {5, 100 <= T <= 1000}}]), I am able to plot the graph nicely. But once I use a straight line equation in that Piecewise function, integration shows weird answers. What am I doing wrong here?

Remove["Global`*"];
f[p_, T_] := 1/(Exp[p/T] + 1)
Needs["Rubi`"]

\[Theta] = (7*10^-11)^(1/2)
Subscript[M, 1] = 2*10^-3
Tf =100

Subscript[C, \[Alpha]][T_] = 
 Piecewise[{{1.2, T < 100}, {8, T > 1000}, {8 + (900/6.8)*(T - 1000), 
    100 <= T <= 1000}}]

\[CapitalGamma][p_, T_] = 
 Subscript[C, \[Alpha]][T]*10^-11 *  p* T^4 // PiecewiseExpand

\[CapitalDelta][p_] := (Subscript[M, 1]^2 /(2*p))

h[p_, T_] = (1/
      8) (\[CapitalGamma][p, T]* \[CapitalDelta][p]^2 * 
       Sin[2*\[Theta]]^2)/(\[CapitalDelta][
         p]^2 + (\[CapitalGamma][p, T]/2)^2) // FullSimplify // 
  PiecewiseExpand

fDW[p_, t_] = 
 f[p*t/Tf, t] Integrate[h[p*t/1000, t], t] // PiecewiseExpand

x[a_, b_] := a/b

xf = Table[{x[p, T], x[p, T]^2*fDW[p, T]}, {p, 1, 10^3}, {T, 100, 500,
     100}] // N

plot1 = ListLogLogPlot[xf, AxesLabel -> {x, Subscript[x\.b2f, DW]}, 
  Ticks -> {Automatic, Automatic}, 
  PlotLabel -> 
   Style["\!\(\*SubscriptBox[\(M\), \(1\)]\)=7keV,\[Theta]\.b2 = \
7x10⁻¹¹"]]






```
$\endgroup$
1
  • $\begingroup$ What "straight line equation" did you use that caused problems? $\endgroup$ – Bob Hanlon Feb 20 at 2:50
2
$\begingroup$

Use exact values until the Table for xf is calculated.

$Version

(* "12.2.0 for Mac OS X x86 (64-bit) (December 12, 2020)" *)

Clear["Global`*"]

f[p_, T_] := 1/(Exp[p/T] + 1)

Use of Rubi is not necessary

(* Needs["Rubi`"] *)

θ = (7*10^-11)^(1/2);
Subscript[M, 1] = 2*10^-3;
Tf = 100;

Subscript[C, α][T_] = Piecewise[
  {{6/5, T < 100}, {8, T > 1000}, {8 + (2250/17)*(T - 1000), 
     100 <= T <= 1000}}];

Γ[p_, T_] = 
  Subscript[C, α][T]*10^-11*p*T^4 // PiecewiseExpand;

Δ[p_] := (Subscript[M, 1]^2/(2*p))

h[p_, T_] = (1/
       8) (Γ[p, T]*Δ[p]^2*
        Sin[2*θ]^2)/(Δ[
          p]^2 + (Γ[p, T]/2)^2) // FullSimplify // 
   PiecewiseExpand;

fDW[p_, t_] = f[p*t/Tf, t] Integrate[h[p*t/1000, t], t] // 
  PiecewiseExpand;

x[a_, b_] := a/b

EDIT: To eliminate complex artifacts (imaginary parts negligible compared to real part) use high precision and Chop[#, 10^-30]&. Also note that the iterators in the Table are done in the opposite order to group by T values .

xf = Table[{x[p, T], x[p, T]^2*fDW[p, T]},
     {T, 100, 500, 100}, {p, 1, 10^3}] //
    N[#, 20] & // Chop[#, 10^-30] &;

plot1 = ListLogLogPlot[xf,
  Joined -> True, AxesLabel -> {x, Subscript[x\.b2f, DW]}, 
  Ticks -> {Automatic, Automatic}, 
  PlotLabel -> 
   Style["\!\(\*SubscriptBox[\(M\), \(1\)]\)=7keV,\[Theta]\.b2 = 7x10⁻¹¹"],
  PlotLegends ->
   LineLegend[Range[100, 500, 100],
    LegendLabel -> "T"]]

enter image description here

$\endgroup$
2
  • $\begingroup$ Thanks for the answer. But even when I am running this, the integration result still contains RootSum and I am getting complex values in the Table for xf where it is supposed to be just real numbers $\endgroup$ – annoying_noob Feb 20 at 10:22
  • 2
    $\begingroup$ There is nothing wrong with RootSum it is just a means of representing an exact number that cannot be represented by radicals. You can eliminate the complex artifacts by using Chop with the threshold set several orders of magnitude below the smallest real value. See edit. $\endgroup$ – Bob Hanlon Feb 20 at 15:26

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.