Differential complex equations in cosmological problem

Question

I want to solve a differential equation of second order with complex initial conditions. Despite I have wrote all the code with detail, I am still facing some problems I might want to review here. The code is as follows:

ode1 = x''[t] + 3 x'[t] Sqrt[(x'[t])^2/6 + (1 - Exp[-x[t] Sqrt[2/3]])^2/4] + Sqrt[3/2] Exp[-x[t] Sqrt[2/3]] (1 - Exp[-x[t] Sqrt[2/3]]) == 0
ode1IC = {x'[0] == -0.008226306418212731, x[0] == 5.630991866033891};
solX = NDSolve[{ode1, ode1IC}, x, {t, 0, 500}][[1, 1]]
ode2 = a'[t]/a[t] == Sqrt[(x'[t])^2/6 + (1 - Exp[-x[t] Sqrt[2/3]])^2/4];
ode2 = ode2 /. solX;
ode2IC = a[0] == 1;
sola = NDSolve[{ode2, ode2IC}, a, {t, 0, 500}][[1, 1]]
ode3 = \[Tau]'[t] == 1/a[t];
ode3 = ode3 /. sola;
ode3IC = \[Tau][149.4517772937791] == 0;
solTao = NDSolve[{ode3, ode3IC}, \[Tau], {t, 0, 500}][[1, 1]]
f = ((7.41193*10^6)^2 - ((3/2 + 1/2*((Sqrt[3/2] Exp[-x[t] Sqrt[2/3]] (1 - Exp[-x[t] Sqrt[2/3]]))/(3/4 (1 - Exp[-x[t] Sqrt[2/3]])^2))^2 + 1/2 ((Exp[-Sqrt[8/3] x[t]] (1 - Exp[Sqrt[2/3] x[t]] (1 - Exp[-Sqrt[2/3] x[t]]))/(3/4 (1 - Exp[-x[t] Sqrt[2/3]])^2))))^2 - 1/4)/(\[Tau][t])^2) /. {solX,solTao};
ode4 = u''[t] + f*u[t] == 0;
ode4 = SetAccuracy[(Solve[ode4, u''[t]] /. Rule -> Equal)[[1, 1]] // 
Simplify, Infinity];
ode4 = ode4 /. {solTao, solX};
ic={u[0]==0.000259728 E^(I (-0.0249668)),u'[0]==1925.09 E^(I (-1.59576))};
ic=SetAccuracy[ic,Infinity];
sol = DSolveValue[{ode4, ic}, u[t], t];

I am interested in Plotting the solution of ode4. When running the Plot for the solution, it shows.

Plot[Abs@sol,{t,0,100}]
DSolveValue::dsvar: The variable specification is not valid.
General::stop: Further output of DSolveValue::dsvar will be suppressed during this calculation.

How do I solve it?

Update:

This is what happens when running the sol code.

NDSolve[{x''[t] + 3 x'[t] Sqrt[(x'[t])^2/6 + (1 - Exp[-x[t] Sqrt[2/3]])^2/4] + Sqrt[3/2] Exp[-x[t] Sqrt[2/3]] (1 - Exp[-x[t] Sqrt[2/3]]) == 0, a'[t]/a[t] == Sqrt[(x'[t])^2/6 + (1 - Exp[-x[t] Sqrt[2/3]])^2/4], \[Tau]'[t] == 1/a[t], u''[t] + ((7.41193*10^6)^2 - ((3/2 + 1/2*((Sqrt[3/2] Exp[-x[t] Sqrt[2/3]] (1 - Exp[-x[t] Sqrt[2/3]]))/(3/4 (1 - Exp[-x[t] Sqrt[2/3]])^2))^2 + 1/2 ((Exp[-Sqrt[8/3] x[t]] (1 - Exp[Sqrt[2/3] x[t]] (1 - Exp[-Sqrt[2/3] x[t]]))/(3/4 (1 - Exp[-x[t] Sqrt[2/3]])^2))))^2 - 1/4)/(\[Tau][t])^2)*u[t] == 0, x'[0] == -0.008226306418212731, x[0] == 5.630991866033891, a[0] == 1, \[Tau][149.4517772937791] == 0, u[0] == 0.000259728 E^(I (-0.0249668)), u'[0] == 1925.09 E^(I (-1.59576))}, {x, \[Tau], a, u}, {t, 0, 100}]
Power::infy: Infinite expression 1/0.^2 encountered.
Infinity::indet: Indeterminate expression 1.5 +ComplexInfinity+ComplexInfinity encountered.
General::stop: Further output of Power::infy will be suppressed during this calculation.
Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.
General::stop: Further output of Infinity::indet will be suppressed during this calculation.
NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.`.

Both methods seem to have problems. Which is the best way to solve them?

Second update:

I figured out that I can solve the system with NDSolve, the trouble is that it takes so much time to solve. I need to make it efficient, but even though with a million steps, I see that it has only been solved from 0 to 0.08, when I want to solve it for the first 100 in $t$.

solode = NDSolve[{u''[t] + u[t]*(k55^2 - f) == 0, u[0] == 0.000259728 E^(I (-0.0249668)), u'[0] == 1925.09 E^(I (-1.59576))}, u, {t, 0, 20}, WorkingPrecision -> 30, MaxSteps -> 1000000]

What is difficult to understand for me is that the first time I ran this code, it achieved to cover a solution from $t=0$ to $t=20$, and now it doesn't go further more from zero. Is there any way to solve this efficiency problem?

Observe here that $k_{55}$ is a fixed constant, while $f$ is a function, which is

\[Nu] = 3/2 + eps1 + 1/2*eps2
f = (\[Nu]^2 - 1/4)/(\[Tau]^2)

The terms that form $f$ are also functions with respect to $t$, all of them well defined and already solved by another NDSolve previously in the code. $u(t)$ is the function that I need a solution for; however, since it oscillates too much in extended times, it can be difficult without precision. Next image is the

ReImPlot[u[t]/.solode,{t,0,20}]

Answer 1

Solving this system of ODEs does not seem feasible, either symbolically or numerically. However, it is feasible to provide approximate solutions for t < 100.

For convenience, simultaneously solve the first three ODEs numerically.

ode1 = x''[t] + 3  x'[t] Sqrt[(x'[t])^2/6 + (1 - Exp[-x[t]  
  Sqrt[2/3]])^2/4] + Sqrt[3/2] Exp[-x[t] Sqrt[2/3]]  (1 - Exp[-x[t]  
  Sqrt[2/3]]) == 0;
ode1IC = {x'[0] == -0.008226306418212731, x[0] == 5.630991866033891};

ode2 = a'[t]/a[t] == Sqrt[(x'[t])^2/6 + (1 - Exp[-x[t] Sqrt[2/3]])^2/4];
ode2IC = a[0] == 1;

ode3 = \[Tau]'[t] == 1/a[t];
ode3IC = \[Tau][149.4517772937791] == 0;

sol123 = NDSolve[{ode1, ode2, ode3, ode1IC, ode2IC, ode3IC}, 
  {x[t], \[Tau][t]}, {t, 0, 500}] // Flatten

Plot[Evaluate[{x[t], \[Tau][t]} /. sol123], {t, 0, 500}, PlotRange -> All, 
  AxesLabel -> {t, "x,\[Tau]"}, LabelStyle -> {10, Bold, Black}]

The small oscillations in x gradually die away with increasing t. a is not plotted, because it is not needed in subsequent calculations. Next, compute and plot f, which appears in ode4.

f = ((7.41193*10^6)^2 - ((3/2 + 1/2*((Sqrt[3/2] Exp[-x[t] Sqrt[2/3]] 
  (1 - Exp[-x[t] Sqrt[2/3]]))/(3/4 (1 - Exp[-x[t] Sqrt[2/3]])^2))^2 + 
  1/2 ((Exp[-Sqrt[8/3] x[t]] (1 -  Exp[Sqrt[2/3] x[t]] 
  (1 - Exp[-Sqrt[2/3] x[t]]))/(3/4 (1 - Exp[-x[t] Sqrt[2/3]])^2))))^2 - 1/4)
  /(\[Tau][t])^2);

Plot[f /. sol123, {t, 0, 145}, AxesLabel -> {t, "x,\[Tau]"}, 
  LabelStyle -> {10, Bold, Black}]

For larger values of t, f is strictly negative and oscillates between -10^17 and -10^30. The analysis below focuses on t < 100. As seen in the plot just above, f is approximately constant for t < 20, and ode4 produces solutions oscillating at a frequency Sqrt[f /. sol123 /. t -> 0], i.e., 7.41193*10^6. It is for this reason that NDSolve takes seemingly forever to solve ode4. On the other hand, ode4 can be solved approximately using the WKB approximation, if desired, for t < 100, except near f = 0. There, an approximate symbolic solution exists:

FindRoot[f /. sol123, {t, 36}][[1, 2]]
(* 36.7625 *)
Series[f /. sol123, {t, %, 1}] // Normal
(* -0.015625 - 7.44662*10^12 (-36.7625 + t) *)

The constant term is neglibile and can be ignored, allowing ode4 to be solved in terms of Airy functions.

as = DSolveValue[u''[t] + %[[2]] u[t] == 0, u[t], t] /. {C[1] -> 1, C[2] -> 0}
(* AiryAi[2.62237*10^-9 (-2.73757*10^14 + 7.44662*10^12 t)] *)

(The second solution, deleted here, grows exponentially for large t.)

Plot[as, {t, 36.7625 - .002, 36.7625 + .002}, AxesLabel -> {t, u}, 
  LabelStyle -> {10, Bold, Black}]

The WKB approximation, which I shall not take the time to apply here in detail, is described in https://en.wikipedia.org/wiki/WKB_approximation. From it, the amplitude of the oscillations varies as f^(-1/4) for t < 36.5. Evaluating f at t = 0 and t = 36.5 indicates that the amplitude of oscillations in u grows by about a factor of

((f /. sol123 /. t -> 36.5)/(f /. sol123 /. t -> 0))^(-1/4)
(* 2.28105 *)