Any ideas on how to use NDSolve to solve this nonlinear 2nd order ODE

Question

I can't get to work this code to solve a nonlinear second order boundary value problem.

ClearAll["Global`*"]; 
a = 1; 
b = 2; 
q1[t_] := {Sin[2*Pi*t], Cos[2*Pi*t]}
q2[t_] := {a*Sin[2*Pi*t], b*Cos[2*Pi*t]}
f[x_, y_] := x . y
eqn = 0.5*Derivative[2][y][t]*f[q2[y[t]], q1[t]] - 
    3*Derivative[1][y][t]^2*
           f[Derivative[1][q2][y[t]], q1[t]] + 
    Derivative[1][y][t]*f[q2[y[t]], Derivative[1][q1][t]] + 
         
    4*Derivative[1][y][t]^2*Sqrt[Derivative[1][y][t]]*
     f[q2[y[t]], q2[y[t]]] == 0; 
bc = {y[0] == 0, y[1] == 1}; 
sol = NDSolve[{eqn, bc}, y, {t, 0, 1}, 
  Method -> "StiffnessSwitching"]
Plot[Evaluate[y[x] /. sol], {x, 0, 1}, PlotRange -> All]

I get immediately the error message

(1) "NDSolve::ndsz: At t == 0.2499999999990644`, step size is effectively zero; singularity or stiff system suspected."

When I add "AccuracyGoal -> 18, PrecisionGoal -> 18" to this code, I get the message

(2) "NDSolve::mxst: Maximum number of 79973 steps reached at the point t == 0.2499927591408232`."

When I reduce to "AccuracyGoal -> 8, PrecisionGoal -> 8", I get message (1).

One difficulty in solving this ODE may be the square root in the equation. So, I converted it into a nonlinear system of first order equations without the square root using the transformation $u=\sqrt{y'}$. Then, $u^2 = y'$ and $2uu'=y''$. Also, let $v=y$, then $v'=u^2$. The 1st order system is: $\{u'-u^3f_1(t,v)+uf_2(t,v)+u^4f_3(t,v)=0,v'=u^2\}$, where $f_1(t,v)=q_2'(y).q_2(t)/q_2(y).q_1(t)$, $f_2(t,v)=q_2(y).q_1'(t)/q_2(y).q_1(t)$, $f_3(t,v)=q_2'(y).q_2(t)/q_2(y).q_1(t)$

The code is:

s = NDSolve[{Derivative[1][u][t] - 
           u[t]^3*(3*(f[Derivative[1][q2][v[t]], q1[t]]/
                   f[q2[v[t]], q1[t]])) + 
           u[t]*(f[q2[v[t]], Derivative[1][q1][t]]/
                f[q2[v[t]], q1[t]]) + u[t]^4*
             (4*(f[Derivative[1][q2][v[t]], q2[v[t]]]/
                   f[q2[v[t]], q1[t]])) == 0, 
       Derivative[1][v][t] - u[t]^2 == 0, v[0] == 0, v[1] == 1}, 
     {u, v}, {t, 0, 1}, Method -> "StiffnessSwitching", 
     AccuracyGoal -> 18, PrecisionGoal -> 18]
Plot[Evaluate[{v[t], u[t]}], {t, 0, 1}, PlotRange -> All]

I get the message, "The scaled boundary value residual error of 4.714045207910324`*^7 indicates that the boundary values are not satisfied to specified tolerances. Returning the best solution found." It returns constant solutions.

Answer 1

Many nonlinear ODEs have no solution and this is such a case, as can be seen from

ff = ParametricNDSolveValue[{eqn, y[0] == 0, y'[0] == y0}, y[1], {t, 0, 1}, {y0}];
Plot[ff[w], {w, 0, 4}, AxesLabel -> {y'[0], y[1]}, 
    LabelStyle -> {12, Bold, Black}, PlotRange -> All]

For no value of y'[0] does y[1] equal 1. Here are some typical solutions.

gg = ParametricNDSolveValue[{eqn, y[0] == 0, y'[0] == y0}, y[t], 
    {t, 0, 1}, {y0}];
Plot[Evaluate@Table[gg[i], {i, {.01, .1, 1, 4}}], {t, 0, 1}, 
    AxesLabel -> {t, y[t]}, LabelStyle -> {12, Bold, Black}, PlotRange -> All]