Using NDSolve and complex differential equations

I am trying to solve this system of equations with NDSolve

 sol = 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 - a''[t]/a[t] + (a[t])^2 (3 (x'[t])^2 - (x'[t])^4/(2*(Sqrt[(x'[t])^2/6 + (1 - Exp[-x[t] Sqrt[2/3]])^2/4])^2) + 2*x'[t]*Sqrt[3/2]*Exp[-Sqrt[2/3]*x[t]] (1 - Exp[-Sqrt[2/3]*x[t]])/Sqrt[(x'[t])^2/6 + (1 - Exp[-x[t] Sqrt[2/3]])^2/4] + 2*Exp[-Sqrt[8/3]*x[t]] - Exp[-Sqrt[2/3]*x[t]]))*u[t] == 0,
u[0] == 0.000259728 E^(I (-0.0249668)),
u'[0] == 1925.09 E^(I (-1.59576)),
\[Tau](149.4517772937791) == 0,
a[0] == 1,
x'[0] == -0.008226306418212731,
x[0] == 5.630991866033891},
{x, a, \[Tau], u}, {t, 0, 500}];


However, it gives me the following error

 NDSolve::ndode: The equations {{-((I E^(-I InterpolatingFunction[{<<1>>},{<<13>>},{<<1>>},{<<3>>},{<<1>>}][t]) <<1>>[t])/Sqrt[2])}[0]==5.63099,({-((I E^<<1>> InterpolatingFunction[{{<<2>>}},<<4>>][t])/Sqrt[2])}^\[Prime])[0]==-0.00822631,Sqrt[3/2] E^(-Sqrt[(2/3)] {Times[<<4>>]}[t]) (1-E^Times[<<3>>])+3 ({-I Power[<<2>>] Power[<<2>>] InterpolatingFunction[<<5>>][<<1>>]}^\[Prime])[t] Sqrt[1/4 Power[<<2>>]+1/6 Power[<<2>>]]+({-((I E^Times[<<2>>] InterpolatingFunction[{<<1>>},{<<13>>},{<<1>>},{<<3>>},{<<1>>}][t])/Sqrt[2])}^\[Prime]\[Prime])[t]==0} are not differential equations or initial conditions in the dependent variables {a,u,\[Tau]}.


Even after trying to modify it, I do not see a solution. It works if I do not add the equations and initial conditions for $$u(t)$$, which is a complex function, since it has complex initial conditions.

At the end, I want to get something like this

The blue function represents the real part while the orange one is the complex branch.

First small syntax error \[Tau](149.4517772937791) == 0 shuld be \[Tau][149.4517772937791] == 0

I noticed first ode only depends on its self x. So this is solved on its own.

Second ode a[t] depends on x. So after solving the x ode, plugin the solution in it and solve for a[t].

Third ode \[Tau] depends on a. So plugin the a[t] solution and solve for \[Tau][t]

The last ode u[t] depends on both x and a. Plugin these solutions in it and now solve for u[t]

Here NDSolve is taking long time for me to wait for it. But it is running and no error. So now you only have one ODE to figure out why NDSolve is taking long time on (it is the one with complex IC.

code

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 = τ'[t] == 1/a[t];
ode3 = ode3 /. sola;
ode3IC = τ[149.4517772937791] == 0;
solTao = NDSolve[{ode3, ode3IC}, τ, {t, 0, 500}][[1, 1]]


ode4 = u''[
t] + ((7.41193*10^6)^2 -
a''[t]/a[
t] + (a[t])^2 (3 (x'[t])^2 - (x'[
t])^4/(2*(Sqrt[(x'[t])^2/
6 + (1 - Exp[-x[t] Sqrt[2/3]])^2/4])^2) +
2*x'[t]*Sqrt[3/2]*
Exp[-Sqrt[2/3]*x[t]] (1 - Exp[-Sqrt[2/3]*x[t]])/
Sqrt[(x'[t])^2/6 + (1 - Exp[-x[t] Sqrt[2/3]])^2/4] +
2*Exp[-Sqrt[8/3]*x[t]] - Exp[-Sqrt[2/3]*x[t]]))*u[t] == 0;
ode4 = ode4 /. {sola, solX};
ode4IC = {u[0] == 0.000259728 E^(I (-0.0249668)), u'[0] ==
1925.09 E^(I (-1.59576))};

solU = NDSolve[{ode4, ode4IC}, u, {t, 0, 500}][[1, 1]]


waiting....may be need to change options and play around with different methods on this one...but no error...

Update

I left the code for u(t) running all night, and it is still running it.

Here is the issue. Your 4th ode has this form. Lets start from the 4th ode, the one with $$u''(t)$$ using the above code.

ode4 = u''[
t] + ((7.41193*10^6)^2 -
a''[t]/a[
t] + (a[t])^2 (3 (x'[t])^2 - (x'[
t])^4/(2*(Sqrt[(x'[t])^2/
6 + (1 - Exp[-x[t] Sqrt[2/3]])^2/4])^2) +
2*x'[t]*Sqrt[3/2]*
Exp[-Sqrt[2/3]*x[t]] (1 - Exp[-Sqrt[2/3]*x[t]])/
Sqrt[(x'[t])^2/6 + (1 - Exp[-x[t] Sqrt[2/3]])^2/4] +
2*Exp[-Sqrt[8/3]*x[t]] - Exp[-Sqrt[2/3]*x[t]]))*u[t] == 0;
ode4 = SetAccuracy[(Solve[ode4, u''[t]] /. Rule -> Equal)[[1, 1]] //
Simplify, Infinity]


When plugging in the solutions found before for $$a(t),x(t)$$ it is really very simple ode. It has the form of simple harmonic oscillator

$$u''(t) + f(t) u(t) = 0$$

With your complex (!) initial conditions.

The problem is that your $$f(t)$$ is very complicated and time dependent. Lets look at your $$f(t)$$. Using earlier code:

ode4 = ode4 /. {sola, solX}
f = ode4[[2]] /. u[t]*any_ :> any; (*extract f(t) *)
Plot[f, {t, 0, 500}]


It looks like it is constant up to $$t=100$$ or so. So lets try to use $$f(t)$$ as constant and solve for very small $$t$$ and see if this helps

So the ode becomes

 u''[t] +5.49367*10^13  u[t] == 0  (*moved things to the left*)


For small $$t$$ since this is approximation of $$f(t)$$. So the ode now

ode4=u''[t]==u[t]*(f/.t->0)
ode4=SetAccuracy[ode4,Infinity]
ic={u[0]==0.000259728 E^(I (-0.0249668)),u'[0]==1925.09 E^(I (-1.59576))};
ic=SetAccuracy[ic,Infinity]


Solve

sol = DSolveValue[{ode4, ic}, u[t], t]


The above is the solution which is valid assuming $$f(t)$$ has the above constant value, so for small $$t$$.

ReImPlot[sol, {t, 0, 0.00001}]


Plot[Abs@sol,{t,0,0.00001}]


NDSolve can now solve it also:

NDSolve[{ode4, ic}, u[t], {t, 0, 1}]


So the main issue I see is that your $$f(t)$$ is very complicated over large time frame. See if you can fit it to simpler function or solve for small $$t$$ only and approximate it to constant.

ps. did your ode come from some physical model? It is little strange to have second order ode with real coefficients having complex initial conditions.

• Do you have any method to solve it quickly? Referring to $u(t)$. I changed $t$-range from 0 to 150, since I am only interested on that interval Dec 21, 2023 at 4:07
• I left the code for $u(t)$ running all night, and it is still running it. I guess it is needed another method Dec 21, 2023 at 13:49
• When I limit the number of steps to 10 million, it shows that NDSolve::mxst: Maximum number of 10000000 steps reached at the point t == 0.20694123245297608, and now I just want to solve $u$ for $t$ from 0 to 50, but it seems to require too much power. Dec 21, 2023 at 17:32
• My model comes from cosmological inflation equations, trying to solve the Mukhanov-Sasaki equation for a Starobinsky potential Dec 22, 2023 at 16:28
• Since I have to get something as the image I have added after an edit, I do not see fixing $f$ as a constant as a solution, since I have to see the behavior around $t=0$ to $t=50$. $\mathcal{R}_k=H(t)/(a(t)ϕ′(t))u_k(t)$, where H[t]=Sqrt[(x'[t])^2/6 + (1 - Exp[-x[t] Sqrt[2/3]])^2/4]. Notice that I took x[t]=\[Phi][t]. My goal is to graph $|\mathcal{R}_k|^2$, but even if I am running NDSolve[{u''[t] + f*u[t] == 0, u[0] == 0.000259728 E^(I (-0.0249668)), u'[0] == 1925.09 E^(I (-1.59576))}, u, {t, 20, 41}], it takes too long to get an answer. Dec 22, 2023 at 19:29