1

我有一个非常复杂的积分,必须进行数值评估。被积函数 (p) 是两个变量 (z, Lx) 的函数,积分在 Lx 上。

我的输出出现错误,从这里我得到的问题是我也需要集成另一个变量(这正是我想要避免的)。

有没有办法解决这个问题?该代码有效,因此如果您尝试运行它,您将看到我想要获得的内容,唯一的问题是结果与文献相比略有错误(因子 3 或 5,取决于所采用的参数)。

这是我的代码,在此先感谢。

constants = {A -> 1.21/10^6, p1 -> 4.62, p2 -> -1.15, zc0 -> 2.49, 
  alpha -> 0.2, La -> 10^45.74, L0 -> 10^44.25, g1 -> 1.01, 
  g2 -> 2.38, \[CapitalOmega]m -> 0.3, \[CapitalOmega]l -> 
   0.7, \[CapitalOmega]r -> 0, h -> 0.7, c -> 300000, H0 -> 70, 
  Fov -> 4.476176190127657*^-7, Mpc2cm -> 3.08567758*10^24}

phizero[Lx_] := A/((Lx/L0)^g1 + (Lx/L0)^g2);
p0[L1_, L2_] := NIntegrate[ Lx*phizero[Lx] /. constants, {Lx, L1, L2}];
zc[Lx_] := 
  Piecewise[{{zc0, Lx >= La}, {zc0*(Lx/La)^alpha, Lx < La}}];
e[z_, Lx_] := 
  Piecewise[{{(1 + z)^p1, 
     z <= zc[Lx]}, {(1 + zc[Lx])^p1*((1 + z)/(1 + zc[Lx]))^p2, 
     z > zc[Lx]}}];
phi[z_, Lx_] := phizero[Lx]*e[z, Lx];
 p[z_, L1_, L2_] := 
  NIntegrate[Lx^-1 phi[z, Lx] /. constants, {Lx, L1, L2}];
(*p[3,10^44.5,10^45.5]*)
Manipulate[
 LogPlot[p[z, 10^L1, 10^L2], {z, 0.1, 7}, GridLines -> Automatic, 
  AxesLabel -> {"z", \[Phi]  [\!\(TraditionalForm\`
\*FractionBox[\(1\), 
SuperscriptBox[\(Mpc\), \(3\)]]\)]} ], {L1, 42, 47}, {L2, 43, 47}]
4

0 回答 0