2
$\begingroup$

I am trying to solve the following differential equations: \begin{align} \dot{a}(\tau) &= \frac{a(\tau)}{\sqrt{3}}\left[ \rho_R(\tau) + \rho_V(\tau) \right]^{1/2}\\[0.25cm] \dot{\rho}_R(\tau) &+ 4\frac{\dot{a}(\tau)}{a(\tau)}\rho_R(\tau)=-\dot{\rho}_V(\tau) \end{align} with the initial conditions $$ a(\tau = 0) = 1,\qquad \rho_R(\tau = 0) = 1. $$ The function $\rho_V(\tau)$ is the source of all of the problems. It is of the form $\rho_V(r) = \exp(-I(\tau))$, where $$ I(\tau) = \frac{4\pi}{3} \int_{0}^{\tau}\Gamma(\tau')a^3(\tau')R^3(\tau', \tau)\ d\tau' $$ and $$ R(\tau', \tau) = \int_{\tau'}^\tau \frac{d\tau''}{a(\tau'')}, \qquad \Gamma(\tau) = \Gamma_0e^{\beta\tau} $$ Let $\Gamma_0 = \beta = 1$ for simplicity.

I am having trouble integrating this code using NDSolve since the integrals depend on $a(\tau)$. I would imagine that this is possible to evaluate since the integration only occurs until the currently evaluated time-step, but I am unsure how to do so. Below is my best attempt, but it doesn't work.

The key issue, is that I am unsure how to pass $a(\tau)$ into the function Ii so that I can perform the integration. Thanks for the help!

Γ[t_] := Γ0*Exp[β*t] ; Γ0 = 1; β = 1;

Rr[𝜏p_?NumericQ, 𝜏_?NumericQ, A_?NumericQ] := NIntegrate[1/A, {t, 𝜏p, 𝜏}];
Ii[𝜏_?NumericQ,  A_] := (4*Pi)/3 NIntegrate[Γ[𝜏p]*A^3 Rr[𝜏p, 𝜏, A]^3, {𝜏p, 𝜏0, 𝜏}];

𝜌V[𝜏_ , A_?NumericQ] = Exp[-Ii[𝜏,  A]];
𝜏0 = 0; 𝜏f = 50;
eqs = {
   a'[𝜏] == a[𝜏]/Sqrt[3]*Sqrt[𝜌R[𝜏] + 𝜌V[𝜏, a[𝜏]]],
   𝜌R'[𝜏] + 4*a'[𝜏]/a[𝜏] 𝜌R[𝜏] == -D[𝜌V[𝜏, a[𝜏]], 𝜏],
   𝜌R[𝜏0] == 1, a[𝜏0] == 1
   };

sols = NDSolve[eqs, {𝜏, 𝜏0, 𝜏f}]; 

Update 1

@Alex Trounev has seemed to produce many solutions, all of which agree with each other. I tried to develop another solution, which seems to get close but isn't quite the same as the three results below.

Much like, @Alex Trounev I defined

\begin{equation} R(\tau) = \int_0^\tau \frac{d \tau'}{a(\tau')}, \qquad \dot{R}(\tau) = \frac{1}{a(\tau)} . \end{equation}

Given $I(\tau)$ above, you can expand the cubic term and then define four variables

\begin{equation} v_i(\tau) = \int_{0}^{\tau}\Gamma(\tau')a^3(\tau) R^i(\tau)\ d\tau \end{equation}

each of which corresponds to the differential equation

\begin{equation} \dot{v}_i(\tau) = \Gamma(\tau')a^3(\tau) R^i(\tau), \qquad i = 0,1,2,3 \end{equation}

This allows us to write the function $I(\tau)$ as

\begin{equation} I(\tau) = \frac{4\pi}{3}\left\{ R^3(\tau)v_0(\tau) - 3R^2(\tau)v_1(\tau) + 3R(\tau)v_2(\tau) + v_3(\tau) \right\} . \end{equation}

The initial conditions for each $v_i(\tau = 0) = 0$. I tried integrating this system (see below). While I get a similar result, the height of the peak is significantly different. It also seems like there is an issue with Gam0 = 1. I'm not particularly sure what is causing the cusp at t ~ 2.

I am not sure where my error is here, so any input would be appreciated.

Gam[t_] := Gam0*Exp[beta*t] ; Gam0 = 1; beta = 1;

t0 = 0; tf = 4;
Ii[t_] := (4 Pi)/
   3 (Rr[t]^3 v0[t]  - 3* Rr[t]^2 v1[t] + 3* Rr[t]*v2[t] - v3[t]);
eqs = {
   a'[t] == a[t]/Sqrt[3]*Sqrt[rhoR[t] + Exp[-Ii[t]]],
   rhoR'[t] + 4*a'[t]/a[t] rhoR[t] == -D[Exp[-Ii[t]], t],
   Rr'[t] == 1/a[t],
   v0'[t] == Gam[t]*a[t]^3,
   v1'[t] == Gam[t]*a[t]^3*Rr[t],
   v2'[t] == Gam[t]*a[t]^3*Rr[t]^2,
   v3'[t] == Gam[t]*a[t]^3*Rr[t]^3,
   rhoR[t0] == 1,
   a[t0] == 1,
   Rr[t0] == 0,
   v0[t0] == 0,
   v1[t0] == 0, 
   v2[t0] == 0,
   v3[t0] == 0 
   };

sols = NDSolve[eqs, {a, rhoR, Rr, v0, v1, v2, v3}, {t, t0, tf},
   WorkingPrecision -> MachinePrecision];

enter image description here

Update 2

My issues was one rouge minus sign in front of v3 (corrected above). The code now reproduces @Alex Trounev's quite well. Problem solved!

enter image description here

$\endgroup$
1

1 Answer 1

3
$\begingroup$

We can solve this problem using method described in my answer here. First we transform model to the system of 3 equations using equation $R'=1/a$. Hence we have \begin{align} \dot{a}(\tau) &= \frac{a(\tau)}{\sqrt{3}}\left[ \rho_R(\tau) + \rho_V(\tau) \right]^{1/2}\\[0.25cm] \dot{\rho}_R(\tau) &+ 4\frac{\dot{a}(\tau)}{a(\tau)}\rho_R(\tau)=-\dot{\rho}_V(\tau) \\[0.25cm] \dot{R}(\tau) &= \frac{1}{a(\tau)} \end{align}
with updated initial conditions

$a(0) = 1,\qquad \rho_R(0) = 1,\qquad R(0)=0.$

We define integral $I(\tau)$ as

$I(\tau) = \frac{4\pi}{3} \int_{0}^{\tau}\Gamma(\tau')a^3(\tau')(R(\tau)-R(\tau'))^3\ d\tau'$

First derivative of this integral is given by

$I_1=I'(\tau) = 4\pi \int_{0}^{\tau}\Gamma(\tau')(a^3(\tau')/a(\tau))(R(\tau)-R(\tau'))^2\ d\tau'$

Then we map $I, I_1$ on $(-1,1)$ using substitution $\tau'=\frac{1}{2} \tau(1+z)$. Last integrals we compute using Gauss quadrature rule, for this we call

 Get["NumericalDifferentialEquationAnalysis`"];

Code to solve the problem for $0\le \tau\le 4$

UE[m_, t_] := EulerE[m, t];
psi[k_, n_, m_, t_] := 
  Piecewise[{{2^(k/2) UE[m, 2^k t - 2 n + 1], (n - 1)/2^(k - 1) <= t <
       n/2^(k - 1)}, {0, True}}];
PsiE[k_, M_, t_] := 
 Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 3; M0 = 6; With[{k = k0, M = M0}, 
 nn = Length[Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]];
dx = 1/(nn); xl = Table[l*dx, {l, 0, nn}]; xcol = 
 Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk = 
 With[{k = k0, M = M0}, PsiE[k, M, t1]]; In1 = 
 With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Psi[x_] := Psijk /. t1 -> x;
in1[x_] := In1 /. t1 -> x; var1 = Array[aa, {nn}]; var2 = 
 Array[b, {nn}]; var3 = Array[c, {nn}];
a[t_] := var1 . in1[t] + a0; a1[t_] := var1 . Psi[t]; 
rhoR[t_] := var2 . in1[t] + b0; rhoR1[t_] := var2 . Psi[t]; 
R[t_] := var3 . in1[t] + c0; R1[t_] := var3 . Psi[t];

np = nn; g = GaussianQuadratureWeights[np, -1, 1]; points = 
 g[[All, 1]];
weights = g[[All, 2]];
Int[ff_, z_] := Sum[(ff /. z -> points[[i]])*weights[[i]], {i, 1, np}];
intNum[t_] := 
  t 2 Pi/3 Int[
    Exp[t/2 (1 + z)] a[t/2 (1 + z)]^3 (R[t] - R[t/2 (1 + z)])^3, z];

int1Num[t_] := 
  t 2 Pi Int[
    Exp[t/2 (1 + z)] (a[t/2 (1 + z)]^3/
       a[t]) (R[t] - R[t/2 (1 + z)])^2, z];

T = 4; eqs = 
 Table[{a1[t]/T == a[t]/Sqrt[3]*Sqrt[rhoR[t] + Exp[-T intNum[t]]], 
    rhoR1[t]/T + 4*a1[t]/a[t] rhoR[t]/T == 
     T int1Num[t] Exp[-T intNum[t]], R1[t]/T == 1/a[t]}, {t, xcol}] //
   Flatten;

Solution

{a0 = 1, b0 = 1, c0 = 0}; var = Join[var1, var2, var3]; sol = 
 FindRoot[eqs, Table[{var[[i]], 1/10}, {i, Length[var]}], 
  MaxIterations -> 10000];

Visualization

{Plot[Evaluate[a[t/T] /. sol], {t, 0, T}, 
  AxesLabel -> {"\[Tau]", "a"}], 
 Plot[Evaluate[rhoR[t/T] /. sol], {t, 0, T}, 
  AxesLabel -> {"\[Tau]", "\[Rho]R"}], 
 Plot[Evaluate[R[t/T] /. sol], {t, 0, T}, 
  AxesLabel -> {"\[Tau]", "R"}]}

Figure 1

Update 1. We can also use method proposed by Ulrich Neumann and adopted for this problem as follows

nl1 = NestList[
   Function[{fa}, 
    Block[{int, int1, fR}, 
     int[tau_?NumericQ] := 
      Block[{s}, 
       4 Pi/3 NIntegrate[
         Exp[s] fa[[1]][s]^3 (fR[s, tau])^3, {s, 0, tau}, 
         Method -> {{Automatic, "LocalAdaptive"}[[-1]], 
           "SymbolicProcessing" -> 0}]]; 
     int1[tau_?NumericQ] := 
      Block[{s}, 
       4 Pi NIntegrate[
         Exp[s] fa[[1]][s]^3/fa[[1]][tau] (fR[s, tau])^2, {s, 0, tau},
          Method -> {{Automatic, "LocalAdaptive"}[[-1]], 
           "SymbolicProcessing" -> 0}]]; 
     fR[s_?NumericQ, tau_?NumericQ] := 
      Block[{s1}, 
       NIntegrate[1/fa[[1]][s1], {s1, s, tau}, 
        Method -> {{Automatic, "LocalAdaptive"}[[-1]], 
          "SymbolicProcessing" -> 0}]];
     NDSolveValue[{ae'[tau] == 
        ae[tau]/Sqrt[3]*Sqrt[rhoRe[tau] + Exp[-int[tau]]], 
       rhoRe'[tau] + 4*ae'[tau]/ae[tau] rhoRe[tau] == 
        int1[tau] Exp[-int[tau]], rhoRe[0] == 1, ae[0] == 1}, {ae, 
       rhoRe}, {tau, 0, 4}]]], {1, 1}, 5]; 

Here we use 5 iterations only since it takes a time. To compare with wavelets method we use

pl={Show[Plot[Evaluate[a[t/T] /. sol], {t, 0, T}, 
  AxesLabel -> {"\[Tau]", "a"}, PlotStyle -> Blue],Plot[Table[nl1[[i]][[1]][t], {i, 4, 5}], {t, 0, T}, PlotRange -> All, 
 PlotStyle -> {Red, Dashed}]], Show[Evaluate[rhoR[t/T] /. sol], {t, 0, T}, 
 AxesLabel -> {"\[Tau]", "\[Rho]R"}, PlotRange -> All], Plot[Table[nl1[[i]][[2]][t], {i, 4, 5}], {t, 0, T}, PlotRange -> All, 
 PlotStyle -> {Red, Dashed}]]}

Figure 2

In this picture we see unexplained discrepancies. Maybe we need third code to compare with.

Update 2. Third code is given by (this is appropriate modification code proposed by flinty)

next[pts_, d_] := 
  With[{np = Normal[pts]}, 
   With[{a = 
      Interpolation[Transpose[{np[[All, 1]], np[[All, 2]]}], 
       InterpolationOrder -> 0], p = pts["Part", -1], 
     rhoR = Interpolation[Transpose[{np[[All, 1]], np[[All, 3]]}], 
       InterpolationOrder -> 0], 
     R = Interpolation[Transpose[{np[[All, 1]], np[[All, 4]]}], 
       InterpolationOrder -> 0]}, 
    With[{x0 = p[[1]], a0 = p[[2]], rhoR0 = p[[3]], R0 = p[[4]]}, 
     With[{int = 
        4 Pi/3 NIntegrate[Exp[s] a[s]^3 (R[x0] - R[s])^3, {s, 0, x0}, 
          Method -> {{Automatic, "LocalAdaptive"}[[-1]], 
            "SymbolicProcessing" -> 0}], 
       int1 = 4 Pi NIntegrate[
          Exp[s] a[s]^3 /a[x0] (R[x0] - R[s])^2, {s, 0, x0}, 
          Method -> {{Automatic, "LocalAdaptive"}[[-1]], 
            "SymbolicProcessing" -> 0}]}, 
      pts["Append", {x0 + d, 
        a0 + d*(a[x0]/Sqrt[3]*Sqrt[rhoR[x0] + Exp[-int]]), 
        rhoR0 + d*(int1 Exp[-int] - 
            4*(a[x0]/Sqrt[3]*Sqrt[rhoR[x0] + Exp[-int]])/a[x0] rhoR[
              x0]), R0 + d/a[x0]}]]]]];


points = CreateDataStructure["DynamicArray"];
points["Append", {0, 1, 1, 0}];
With[{n = 200}, result = Nest[next[#, 4./n] &, points, n]];

sol1 = Normal[result]; 

As we can see next is the Euler step implementation. Visualization

pl1 = ListPlot[Transpose[{sol1[[All, 1]], sol1[[All, 2]]}], 
  Frame -> True, PlotStyle -> {Green, PointSize[.005]}];
pl2 = ListPlot[Transpose[{sol1[[All, 1]], sol1[[All, 2]]}], 
  Frame -> True, PlotStyle -> {Green, PointSize[.005]}];
{Show[pl[[1]],pl1],Show[pl[[2]],pl2]}

Figure 3 The resume is that wavelets method probably has large error since we use 24 collocation points only.

With the last code we can compute solution up to $\tau =50$ as follows

With[{n = 500}, result = Nest[next[#, 50./n] &, points, n]];

Visualization Figure 4

$\endgroup$
0

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.