Unable to solve a differential equation system

Question

I fail to solve a system of differential equations with Mathematica. I am quite sure that the system itself is solvable, since I successfully got the result when all variables were dependent on time 't'. But after I change the variable from 't' to 'Ne' (where Ne:=\int H dt) and try to solve it again, Mathematica starts to give me error message. It seems that the 'H[Ne]^2' in eq2 and eq3 is causing the problem, but I cannot further simplify my equations.

I'd appreciate it for any suggestions or comments about what I should do in order to improve my code.

a = 4.10514*10^-6;
W[Ne_] = a/4*y[Ne]^4*Exp[-x[Ne]];
eq1[Ne] = Simplify[H[Ne]^2 - 1/3*W[Ne]];
eq2[Ne] = Simplify[x''[Ne]*H[Ne]^2 + 3*x'[Ne]*H[Ne]^2 + D[W[Ne], x[Ne]]];
eq3[Ne] = Simplify[y''[Ne]*H[Ne]^2 + 3*y'[Ne]*H[Ne]^2 + Exp[x[Ne]]*D[W[Ne], y[Ne]]];
sol1 = NDSolve[{ eq1[Ne] == 0, eq2[Ne] == 0, eq3[Ne] == 0, x[0] == 4, y[0] == 0.56, H[0] == Sqrt[1/3*W[0]], x'[0] == 0, y'[0] == 0}, {H, x, y}, {Ne, 0, 30}, MaxSteps -> 10^8, AccuracyGoal -> 18]

Answer 1

Differentiating eq1 and changing nothing else yields a complete solution:

eq1[Ne] = D[Simplify[H[Ne]^2 - 1/3*W[Ne]], Ne];

Scaled plots:

With[{funcs = {1, 10, 10^5} Through[{x, y, H}[t]]}, 
 Plot[funcs /. First@sol1 // Evaluate, {t, 0, 3},
  PlotLegends -> funcs]
 ]

What difference does it make? If eq1 is an algebraic equation, then the system is a DAE and there only only one numerical integrator that NDSolve can use, namely "IDA". Also the algebraic equation H[Ne]^2 - 1/3*W[Ne] == 0 must be satisfied. There seems to be a problem getting started, the mysterious error tests failure NDSolve::nderr. If we remove AccuracyGoal -> 18, which is probably unachievably high for this system, then integration starts but runs into an NDSolve::ndsz error at Ne == 0.02772779893157676. As Ne approaches 0.02772779893157676, H[Ne] approaches 0. Since H[Ne] is determined by H[Ne]^2..., two things happen: The positive and negative branches of the square root approach each other (which will the numerical solver choose, when they are very close?), and H[Ne] begins to change rapidly ($\sqrt{u}$ has a vertical tangent at $u=0$). NDSolve halts because it can't keep the estimated error low.

By differentiating eq1, H[Ne] is now determined by its derivative. The derivative still goes to infinity when H[Ne] goes to zero, but somehow the stiff solver in LSODA handles it, with what it estimates is an acceptable amount of error. Note that mathematically, it is ambiguous what should happen at H[Ne] == 0. Both the positive and negative square roots are valid solutions. In this case, it so happens that the numerical solver goes for the positive H[Ne], which may be pleasing if H[Ne] represents a physical magnitude. One thing I did not investigate is whether H[Ne] was actually approaching zero as a limit, or it just got close enough that the numerical solver "thought" it was. The computed values only get down to around 10^-8 (or 10^-11 if WorkingPrecision -> 24). If it is supposed to cross zero every oscillation, then there is work to be done.

Answer 2

The problem (besides the non-ideal syntax/style of the code) is an invalid initial condition/ill conditioned system. There is no differential equation for H[Ne] just an algebraic one (eq1[Ne]) it is not necessary to specify a trivial initial condition for it and running H[Ne] as a function in NDSolve is not necessary and might cause problems. What one should do is solve eq1[Ne]==0 for H[Ne]^2 (which is simply H[Ne]^2 -> 1/3*W[Ne]) and substitute the result in eq2[Ne] and eq3[Ne] to eliminate H[Ne]^2 completely. If H[Ne] is a desired quantity on the ODE solution it can be computed from x[Ne] and y[Ne] using W[Ne]. I kept the variable names and style of the original question (while not recommending to code this way in Mathematica):

W[Ne_]=a/4*y[Ne]^4*Exp[-x[Ne]];
eq2[Ne]=Simplify[x''[Ne]*H[Ne]^2+3*x'[Ne]*H[Ne]^2+D[W[Ne],x[Ne]]];
eq3[Ne]=Simplify[y''[Ne]*H[Ne]^2+3*y'[Ne]*H[Ne]^2+Exp[x[Ne]]*D[W[Ne],y[Ne]]];
odes={eq2[Ne]==0,eq3[Ne]==0,x[0]==4,y[0]==0.56,x'[0]==0,y'[0]==0}/.H[Ne]^2->1/3*W[Ne]//Simplify

results in the ODE system:

Note that both ODEs are proportional to y[Ne] and thus trivially true for y[Ne]=0 allowing in principle for arbitrary x[Ne] and derivatives at/beyond this point. This is also the reason why the system is only well defined (just based on the given equations without knowing what they are supposed to describe model) for y[Ne]!=0. The system becomes unstable when y[Ne] approaches zero and hence I modified NDSolve to stop the integration before zero is reached. If one omits this step numerical integration breaks down at this point. The following code

a=4.10514*10^-6;
NeFinal=30;
NDSolve[odes~Join~{WhenEvent[Abs[y[Ne]]<10^-8,NeFinal=Ne;"StopIntegration"]},{x,y},{Ne,0,NeFinal},MaxSteps->10^8,AccuracyGoal->18,Method->"BDF"]//First
{x[Ne],y[Ne]}/.%;
Plot[%,{Ne,0,NeFinal}]
ClearAll[a,NeFinal]

produces the solution: