# Numerically integrating of solution in NDSolve

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*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*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]; 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!

## 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*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"}]} 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[][s]^3 (fR[s, tau])^3, {s, 0, tau},
"SymbolicProcessing" -> 0}]];
int1[tau_?NumericQ] :=
Block[{s},
4 Pi NIntegrate[
Exp[s] fa[][s]^3/fa[][tau] (fR[s, tau])^2, {s, 0, tau},
"SymbolicProcessing" -> 0}]];
fR[s_?NumericQ, tau_?NumericQ] :=
Block[{s1},
NIntegrate[1/fa[][s1], {s1, s, tau},
"SymbolicProcessing" -> 0}]];
NDSolveValue[{ae'[tau] ==
ae[tau]/Sqrt*Sqrt[rhoRe[tau] + Exp[-int[tau]]],
rhoRe'[tau] + 4*ae'[tau]/ae[tau] rhoRe[tau] ==
int1[tau] Exp[-int[tau]], rhoRe == 1, ae == 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]][][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]][][t], {i, 4, 5}], {t, 0, T}, PlotRange -> All,
PlotStyle -> {Red, Dashed}]]} 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[], a0 = p[], rhoR0 = p[], R0 = p[]},
With[{int =
4 Pi/3 NIntegrate[Exp[s] a[s]^3 (R[x0] - R[s])^3, {s, 0, x0},
"SymbolicProcessing" -> 0}],
int1 = 4 Pi NIntegrate[
Exp[s] a[s]^3 /a[x0] (R[x0] - R[s])^2, {s, 0, x0},
"SymbolicProcessing" -> 0}]},
pts["Append", {x0 + d,
a0 + d*(a[x0]/Sqrt*Sqrt[rhoR[x0] + Exp[-int]]),
rhoR0 + d*(int1 Exp[-int] -
4*(a[x0]/Sqrt*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[],pl1],Show[pl[],pl2]} 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]];
`