NDSolve stiff system of ODEs

Question

I am having a hard time trying to solve the following system of ODEs

X'' = e^{2X} + e^{2Y}
Y'' = e^{2Y} - 3 e^{2X}

with initial conditions that mimic a linear behavior of the solutions at $-\infty$. I am using for instance something like

X(u0) = lambdaX * u0
Y(u0) = lambdaY * u0
X'(u0) = lambdaX
Y'(u0) = lambdaY

with $u0 = -10$. The actual values of the parameters $\lambda_X$ and $\lambda_Y$ do not matter much, as long as the solutions are linear at $-\infty$ (ie. both the exponentials have to go to $0$ at $-\infty$, so that $X'' \rightarrow 0$ and $Y'' \rightarrow 0$). I tried for instance with $\lambda_Y = 1.5$ and $\lambda_X = 0.5$.

The basic code that I'm using is something like:

lambdaY = 1.5;
lambdaX = 0.5;
u0 = -10;
sol = NDSolve[{X''[u] == Exp[2 Y[u]] + Exp[2 X[u]], Y''[u] == Exp[2 Y[u]] - 3 Exp[2 X[u]], Y[u0] == u0*lambdaY, X[u0] == u0*lambdaX, Y'[u0] == lambdaY, X'[u0] == lambdaX}, {Y, X}, {u, u0, 10}, Method -> "StiffnessSwitching", AccuracyGoal -> \[Infinity], MaxSteps -> 5000000, InterpolationOrder -> All]

The problem is the usual step size is effectively zero; singularity or stiff system suspected.: the solutions are diverging suddenly at some point and I could not understand what is happening there. I have tried with the standard ways to fix the issue, ie. playing with WorkingPrecision, AccuracyGoal and PrecisionGoal and changing the Method to StiffnessSwitching but none of these worked. I also tried to redefine the variables like

R = e^{2X}
S = e^{2Y}

but I did not manage to avoid the problem.

I am looking for a solution that is asymptotically linear at $-\infty$ for both functions $X$ and $Y$. I am sure that the solution is a simple one and I am missing something. Any suggestion on how to solve this problem?

Thanks a lot.

Answer 1

If we first consider $$X''(u)=e^{2 X(u)} \,,$$ which has a solution $$\log \left(\frac{1}{2} c_1 \left| \sec \left(\frac{1}{2} c_1 (u-{c_2})\right)\right| \right) \,,$$ we see it approaches infinity in finite time.

In the OP's system, $$X''(u)=e^{2 X(u)}+e^{2 Y(u)} > e^{2 X(u)} \,.$$ We should expect that once $X(u)$ starts increasing, it shortly goes to infinity, too.

This can be seen in the plots of $X$, $Y$ and their derivatives:

GraphicsGrid@Map[
  ListLinePlot[# /. First@sol, PlotLabel -> #[u], Frame -> True] &,
  {{X, Y}, {X', Y'}, {X'', Y''}}, {2}]

The orders of growth of $X$, $X'$, $X''$ near the singularity are consistent with a logarithmic pole, too, similar to the solution of the ODE example that began this answer.

Here's a way to see the orders of growth at the singularities $u=u_2 \approx -0.0018$:

With[{u2 = Last@First[X@"Domain" /. First@sol]},
 GraphicsGrid@MapThread[
   Function[{xy, scale},
    Plot[scale #[u] /. First@sol, {u, -0.05, u2}, PlotLabel -> #[u], 
       Frame -> True] & /@ xy],
   {{{X, Y}, {X', Y'}, {X'', Y''}},
    {1/Log[u2 - u], u2 - u, (u2 - u)^2}}]
 ]

$X$, $X'$, $X''$ are asymptotic to $-\log(u_2-u)$, $1/(u_2-u) $, $1/(u_2-u)^2$ respectively; $Y$, $Y'$, $Y''$ are asymptotic to $3\log (u_2-u) $, $-3/(u_2-u) $, $-3/(u_2-u)^2$ respectively.