1
$\begingroup$

Consider an integral

$$ I[f_{n},T]=\int \limits_{0}^{\infty} dE_{n}F(f_{n}(E_{n}),E_{n},T)+\int\limits_{0}^{\infty} dE_{n}\int_{\mathcal{F}(E_{n})}dE_{n}^{'}G(f_{n}(E_{n}),f_{n}(E_{n}^{'}),E_{n},E_{n}^{'},T), $$ where $T>0$ is some parameter, while $f_{n}$ is the distribution function. The integral has the property that for $f_{n} = e^{-E_{n}/T_{n}}$, it 1) is negative if $T_{n}>T$, 2) vanishes at $T_{n} = T$, 3) positive if $T<T_{n}$.

The integral $I$ enters some partial integrodifferential equation on $f_{n}(E_{n},t)$ which I want to solve as a set of ODEs with the help of the discretization of $E$ domain (see this question). Therefore, I need to discretize this integral.

The integral in continuous and discrete forms, $f_{n} = e^{-E_{n}/T_{n}}$, is given below (the final functions are CollisionIntegralIntegratedBoltzmannContinuous[T, Tn], CollisionIntegralIntegratedBoltzmannDiscretized[T, Tn]). To discretize it, I have just used $E_{n} \to E_{\text{min}}+i\Delta E$, where $0\leqslant i \leqslant imax$, and $\int dE_{n} \to \Delta E \sum_{i}$.

(*Some parameters*)
hbar = 6.582119*10^-25;
sToGeVminusOne = 1/hbar;
sToMeVminusOne = sToGeVminusOne*10^-3;
(*Discretization of En values*)
imax = 100;
DEvalue = 0.25;
Emin = 0.01;
Estep[k_] = Emin + DEvalue*k;
fnBoltzmannDiscrete[k_, Tn_] = Exp[-Estep[k]/Tn];
fnThermalBoltzmann[E\[Nu]_, Tn_] = Exp[-E\[Nu]/Tn];


(*Integral I-continuous*)
dS11dE2[E1_, E2_] = -0.0034976546222801287` E1 E2^3*fn[E1] fn[E2];
S12[E1_, T_] = 0.02098592773368077` Exp[-E1/T] E1 T^4;
S21[E1_, T_] = -0.08394371093472308` E1 T^4 fn[E1];
dS221dE2[E1_, E2_, T_] = 
  0.03147889160052116`/(E1^2) Exp[-E1/
     T] T^2 (-E1^2 (E2^2 + 2 E2 T - 2 (-1 + Exp[E2/T]) T^2) - 
     2 E1 T (E2^2 + 2 E2 T - 2 (-1 + Exp[E2/T]) T^2) + 
     2 T^2 ((-1 + Exp[E2/T]) E2^2 - 2 (1 + Exp[E2/T]) E2 T + 
        4 (-1 + Exp[E2/T]) T^2)) fn[E2];
dS222dE2[E1_, E2_, T_] = 
  0.03147889160052116`/(E1^2) Exp[-E1/
     T] T^2 (2 (-1 + Exp[E1/T]) T^2 (E2^2 + 2 E2 T + 4 T^2) - 
     E1^2 (E2^2 + 2 E2 T - 2 (-1 + Exp[E1/T]) T^2) - 
     2 E1 T (E2^2 + 2 E2 T + 2 (1 + Exp[E1/T]) T^2)) fn[E2];
CollisionIntegralContinuous[E1_, T_] := 
 S12[E1, T] + S21[E1, T] + 
  NIntegrate[dS11dE2[E1, E2], {E2, 0, Infinity}] + 
  NIntegrate[dS221dE2[E1, E2, T], {E2, 0, E1}] + 
  NIntegrate[dS222dE2[E1, E2, T], {E2, E1, Infinity}]
CollisionIntegralIntegratedBoltzmannContinuous[T_, Tn_] := 
 NIntegrate[((S12[E1, T] + S21[E1, T]) /. {fn[E1] -> 
        fnThermalBoltzmann[E1, Tn]}) E1^3/(2*Pi^2), {E1, 0, 
    Infinity}] + 
  NIntegrate[(dS11dE2[E1, 
       E2] /. {fn[E1] -> fnThermalBoltzmann[E1, Tn], 
       fn[E2] -> fnThermalBoltzmann[E2, Tn]}) E1^3/(2*Pi^2), {E1, 0, 
    Infinity}, {E2, 0, Infinity}] + 
  NIntegrate[(dS221dE2[E1, E2, 
       T] /. {fn[E1] -> fnThermalBoltzmann[E1, Tn], 
       fn[E2] -> fnThermalBoltzmann[E2, Tn]}) E1^3/(2*Pi^2), {E1, 0, 
    Infinity}, {E2, 0, E1}] + 
  NIntegrate[(dS222dE2[E1, E2, 
       T] /. {fn[E1] -> fnThermalBoltzmann[E1, Tn], 
       fn[E2] -> fnThermalBoltzmann[E2, Tn]}) E1^3/(2*Pi^2), {E1, 0, 
    Infinity}, {E2, E1, Infinity}]

(*Integral I-discretized*)
dS11dE2discrete[i_, j_] = 
  dS11dE2[E1, E2] /. {fn[E1] -> fn[i][t], 
     fn[E2] -> fn[j][t]} /. {E1 -> Estep[i], E2 -> Estep[j]};
S12discrete[i_] = 
  S12[E1, T] /. {fn[E1] -> fn[i][t], fn[E2] -> fn[j][t]} /. {E1 -> 
     Estep[i], E2 -> Estep[j]};
S21discrete[i_] = 
  S21[E1, T] /. {fn[E1] -> fn[i][t], fn[E2] -> fn[j][t]} /. {E1 -> 
     Estep[i], E2 -> Estep[j]};
dS221dE2discrete[i_, j_] = 
  dS221dE2[E1, E2, T] /. {fn[E1] -> fn[i][t], 
     fn[E2] -> fn[j][t]} /. {E1 -> Estep[i], E2 -> Estep[j]};
dS222dE2discrete[i_, j_] = 
  dS222dE2[E1, E2, T] /. {fn[E1] -> fn[i][t], 
     fn[E2] -> fn[j][t]} /. {E1 -> Estep[i], E2 -> Estep[j]};
CollisionIntegralDiscrete[
  i_] := (S12discrete[i] + S21discrete[i] + 
    DEvalue*Sum[dS11dE2discrete[i, j], {j, 0, imax, 1}] + 
    DEvalue*Sum[dS221dE2discrete[i, j], {j, 0, i, 1}] + 
    If[i == imax, 0, 
     DEvalue*Sum[
       dS222dE2discrete[i, j], {j, i + 1, imax, 1}]]) /. {T -> Tg[t]}
CollisionIntegralIntegratedBoltzmannDiscrete[Tg_, 
  Tn_] := (DEvalue*
     Sum[((S12discrete[i] + S21discrete[i]) /. {fn[i][t] -> 
           fnBoltzmannDiscrete[i, Tn]}) (Estep[i])^3/(2*Pi^2), {i, 0, 
       imax, 1}] + 
    DEvalue^2*
     Sum[(dS11dE2discrete[i, 
          j] /. {fn[i][t] -> fnBoltzmannDiscrete[i, Tn], 
          fn[j][t] -> fnBoltzmannDiscrete[j, Tn]}) (Estep[i])^3/(
       2*Pi^2), {i, 0, imax, 1}, {j, 0, imax, 1}] + 
    DEvalue^2*
     Sum[(dS221dE2discrete[i, 
          j] /. {fn[i][t] -> fnBoltzmannDiscrete[i, Tn], 
          fn[j][t] -> fnBoltzmannDiscrete[j, Tn]}) (Estep[i])^3/(
       2*Pi^2), {i, 0, imax, 1}, {j, 0, i, 1}] + 
    DEvalue^2*
     Sum[If[i == imax, 0, 
       Sum[(dS222dE2discrete[i, 
            j] /. {fn[i][t] -> fnBoltzmannDiscrete[i, Tn], 
            fn[j][t] -> fnBoltzmannDiscrete[j, Tn]}) (Estep[i])^3/(
         2*Pi^2), {j, i + 1, imax, 1}]], {i, 0, imax, 1}]) /. {T -> Tg}

Comparing the values of CollisionIntegralIntegratedBoltzmannContinuous[T, Tn], CollisionIntegralIntegratedBoltzmannDiscretized[T, Tn] at $T = T_{n} = 2$, $T = T_{n} = 3$ for DEvalue = 0.25 and imax = 100, I find

Timing[CollisionIntegralIntegratedBoltzmannDiscrete[3, 3]]
CollisionIntegralIntegratedBoltzmannContinuous[3, 3]

{0.941387,-43.7787}

-9.28833*10^-6

This is a very poor agreement. It may be improved by increasing imax and decreasing DEvalue. Namely, for DEvalue = 0.05 and imax = 1000, I get

{98.1421, -0.182558}

The agreement is somewhat better now. However, this improvement is accompanied by a significant increase in timing. This problem suggests that the way of the discretization that is used in my code is poor.

Could you please tell me how to discretize more efficiently, assuming that it may be used for "arbitrary" $f_{n}$ (in the sense that this integral $I$ is evaluated for dynamically changed $f_{n}$)?

$\endgroup$

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.