# Solving Integrate of InterpolatingFunction from an NDSolve

I am interested in solving a system of Differential Equations, particularly

u''[\[Tau]]=-f[t]u[\[Tau]]


This system has for initial conditions $$u(0)=1/\sqrt{2\cdot k_{55}}$$, $$u'(0)=-i\cdot k_{55}/\sqrt{2\cdot k_{55}}$$. As you may have noticed, $$f$$ is in terms of $$t$$, and that's because $$\tau$$ is a function in terms of $$t$$.

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],
τ'[t] == 1/a[t],
x'[0] == -0.008226306418212731,
x[0] == 5.630991866033891,
a[0] == 1,
τ[149.4517772937791] == 0},
{x, τ, a},
{t, 0, 500}]
a = a[t] /. sol;
\[Tau] = \[Tau][t] /. sol;
x = x[t] /. sol;
xt = x'[t] /. sol;
att = a''[t] /. sol;
tEnd = t /. FindRoot[(att == 0), {t, 0, 150}]
aEnd = a /. t -> tEnd;
H = Sqrt[(xt)^2/6 + (1 - Exp[-x*Sqrt[2/3]])^2/4];
V = 3/4 (1 - Exp[-Sqrt[2/3] x])^2;
Vt = Sqrt[3/2] Exp[-x Sqrt[2/3]] (1 - Exp[-x Sqrt[2/3]]);
Vtt = -Exp[-2 Sqrt[2/3]*x]*(Exp[Sqrt[2/3]*x] - 2);
t55 = t /. FindRoot[(aEnd/a == E^55), {t, 0, 150}];
\[Tau]55 = \[Tau] /. t -> t55;
a55 = a /. t -> t55;
eps1 = (1/2) (Vt/V)^2;
eta = Vtt/V;
eps2 = -4 eps1 + 2 eta;
k55 = a55*(H /. t -> t55);
\[Nu] = 3/2 + eps1 + 1/2*eps2;
f = k55^2 - (\[Nu]^2 - 1/4)/(\[Tau]^2);
solu = NDSolve[{u''[t] == -u[t]*f, u[t55] == Exp[I*k55*\[Tau]55]/Sqrt[2*k55], u'[t55]
== -I*k55*Exp[I*k55*\[Tau]55]/(a55*Sqrt[2*k55])}, u[t], {t, t55, tEnd}]


My problem here is that, when trying to get solu (run the NDSolve for $$u(t)$$), it takes too much time to even give a solution. Next, I show the plot of $$f$$ from $$t=0$$ to $$t=t_{\text{End}}$$

Anyway, I am not interested in a solution on that interval but from $$t=t_{55}$$ to $$t=t_{\text{End}}$$, so the initial conditions would be u[t]=Exp[-I*k55*\[Tau]55]/Sqrt[2*k55], and u'(t_{55})=-I*k55*Exp[-I*k55*\[Tau]55]/Sqrt[2*k55]. I hope you can recommend any method to get a solution out of the proposed differential equation, since I need the plot of it. The solution has to look like this:

It only gets to 22 using

solu = NDSolve[{u''[t] == -u[t]*f, u[0] == 1/Sqrt[2*k55],
u'[0] == -I*k55/(a55*Sqrt[2*k55])}, u[t], {t, 0, 22},
Method -> "StiffnessSwitching"]

NDSolve::mxst: Maximum number of 10000 steps reached at the point t ==
22.26269164135809.


Using

solu = NDSolve[{u''[t] == -u[t]*f, u[t55] == Exp[-I*k55*\[Tau]55]/Sqrt[2*k55], u'[t55] == -I*k55*Exp[-I*k55*\[Tau]55]/(a55*Sqrt[2*k55])}, u[t], {t, t55, tEnd}]


Gives the following error: NDSolve::ndsz: At t == 33.42969888667501', step size is effectively zero; singularity or stiff system suspected. And when using MaxSteps -> 10^6 it shows a solution from $$t=t_{55}$$ to $$t = 33.429738218477798903'20$$.

Update:

I took the anzats $$u(t)=e^{\int v(t)dt}$$, therefore, I got that $$v'+v^2+f=0$$. I got $$v$$ using

solv = NDSolve[{v'[t] + (v[t])^2 + f == 0, v[0] == -I*k55}, v[t], {t, 0, 150}]


However, when doing ur = Integrate[v, t], where v=v[t]/.solv, it does not work and does not give any integral at all, not even ur = NIntegrate[v, {t,0,x}]. How do I solve it to get it? With this I will finally get $$u(t)$$.

• It's going to take several steps per oscillation. Looks like lots of oscillations. So maybe MaxSteps -> 10^6 or higher. Commented Jan 27 at 3:46
• @Goofy Yeah, but the think is that I already tried with MaxSteps -> 10^6, and even with MaxSteps -> 10^9. However, I think it is too much for my laptop, and Mathematica just closes randomly. That's why I would like to know if there is any other method to simplify it. Commented Jan 27 at 15:37
• There have been several problems like this (highly oscillatory) posted here, but I don't recall any being solved. You would think such a common problem would have been researched. Maybe try on one of the math sites, scicomp or math or mathoverflow. If there is a math solution, maybe it can be implemented in Mathematica. Commented Jan 27 at 16:34
• f as computed in your question, is a complicated expression enclosed in curly brackets, in other words it is a List. As a result, using it in NDSolve produces unpredictable results. Use First[f] to correct this problem. With that correction , you will see that NDSolve[{u''[t] == -u[t]*First[f], u[0] == 1/Sqrt[2*k55], u'[0] == -I*k55/(a55*Sqrt[2*k55])}, u[t], {t, 0, 10^-5}] // Flatten produces a rapidly oscillating solution of approximately uniform amplitude near t = 0`, as it should. Commented Jan 29 at 14:24
• If you hope to have readers help you, you should first correct known errors in your code. By the way, you question largely duplicates mathematica.stackexchange.com/q/296808/1063. Commented Jan 29 at 15:23