# Differential complex equations in cosmological problem

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}]


• I want to solve a differential equation of second order with complex initial conditions. it will be easier if you first show the actual ODE with its IC you want to solve. Then show your method you tried to use (which you did). But you never showed what is the original ODE with its IC you actually want to solve. Is it ode1 in your code that you want to solve? If so, its IC's are not complex. This is confusing. Also you have typo in pasting your code for f it has -- there near the end which gives syntax error. Commented Jan 22 at 7:28
• You have a double minus character in "f=...". Further, it looks like DSolveValue is unable to solve the DE. Commented Jan 22 at 8:57
• It looks like you're trying to symbolically integrate (DSolveValue) a numerical solution. I thought that usually would not work. Does DSolveValue actually give you a solution? Commented Jan 24 at 16:16
• f is of order 5 10^13 for small t, so u oscillates very rapidly there, which is why the numerical calculation of u is so slow. Certainly, this is inconsistent with the plot in the question. Commented Jan 24 at 17:09
• @bbgodfrey do you think solving the imaginary and real parts of u[t] would make it easier for Mathematica to solve the system? Commented Jan 24 at 17:23

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 *)

• Nice solution to this. However, I do not understand where can I use the proposed initial conditions for the ode4. Do I have to write them in the C[1] and C[2]? Commented Jan 25 at 1:20
• You could use the ode4 boundary conditions to initialized the WKB approximation, then match that solution to the Airy solution near f = 0. However, with of order a million oscillations over that range, the oscillation phase probably would be quite uncertain by time f = 0 is reached. What are you actually trying to determine, the form of the solution or its precise value at some point in t? Commented Jan 25 at 1:42
• I am trying to match the behavior of this solution with another one that I already know. It has to be similar to the last plot in my post, but using the initial conditions I proposed, I got -(-1)^(2/3) 6^(1/3)WeierstrassP[(-(1/6))^(1/3) ((0.000259647 - 6.4839*10^-6 I) + t), {0, -48.0523 - 1924.49 I}] as a solution, but its plot is not what I expected, since the perturbations should follow a Sine, Cosine behavior that is decreasing as t grows. Commented Jan 25 at 1:53
• @JulianYussef The shape of the curve will be sinusoidal, with the amplitude slowly growing until f = 0 is reached, then it will drop exponentially toward zero as shown in my last plot. (Conceivably, it also could grow exponentially from there, but that does not sound physically realistic.) Qualitatively, the solution depends on the shape of f and will differ significantly from my prediction only if f increases with t instead of passing through zero. Commented Jan 25 at 2:01
• Yeah, I got your exact solution, however when I put the initial conditions, it seems that the solution goes crazy with the WeierstrassP function, which is something unusual since I would expect the same sinusoidal behavior as the one you described and posted. DSolveValue[u''[t] + %[[2]] u[t] == 0, u[t], t] /. {C[1] -> 0.000259728 E^(I (-0.0249668)), C[2] -> 1925.09 E^(I (-1.59576))}` Commented Jan 25 at 2:05