Problem with C[i] of a second-order nonlinear ordinary differential equation

Question

I'm trying to solve this differential equation (depending on the K0 parameter):

F[z_] := Sinh[z];
DSolve[{Q''[t] + Sinh[Q[t]] == 0, Q'[0] == 0, Q[0] == Exp[K0 I Pi]}, Q, t];

with no positive results...

DSolve::bvfail: For some branches of the general solution, unable to solve the conditions.
DSolve::bvfail: For some branches of the general solution, unable to solve the conditions.

Now, if I remove the initial conditions, everything works fine, and I get the double result

{{Q -> Function[{t}, -2 I JacobiAmplitude[1/2 Sqrt[-(-2 + C[1]) (t + C[2])^2],
                                          -(4/(-2 + C[1]))]]},
 {Q -> Function[{t}, 2 I JacobiAmplitude[1/2 Sqrt[-(-2 + C[1]) (t + C[2])^2],
                                         -(4/(-2 + C[1]))]]}}

How can I proceed in order to get the two coefficient C[1], C[2] determined with the two initial conditions

Q'[0] == 0, Q[0] == Exp[K0 I Pi]

I tried with Solve and Reduce:

Q[t_] = Q[t] /. First @ EY
Out[507]: -2 I JacobiAmplitude[1/2 Sqrt[-(-2 + C[1]) (t + C[2])^2], -(4/(-2 + C[1]))]
Reduce[Q[0] == E^(I K0 Pi), C[1]]
   Reduce::nsmet: This system cannot be solved with the methods available to Reduce.
Out[508]:Reduce[-2 I JacobiAmplitude[1/2 Sqrt[-(-2 + C[1]) C[2]^2], -(4/(-2 + C[1]))] ==
                E^(I K0 Pi), C[1]]

but they're not working as well... Any ideas?

Answer 1

Once again, Mathematica is not very good at handling elliptic integrals and elliptic functions.

Here is the solution sought, starting from Mathematica's generic solution given in the OP, and helping Mathematica a fair bit in inverting what it couldn't invert:

qfun[t_, K0_] := 2 I JacobiAmplitude[I Sinh[Exp[I K0 π]/2]
                                     (t - EllipticK[-Sinh[Exp[I K0 π]/2]^2]),
                                     -Csch[Exp[I K0 π]/2]^2]

This corresponds to setting C[2] -> -EllipticK[(2 - C[1])/4] and C[1] -> 2 Cosh[Exp[I K0 π]] in the generic solution.

Plot it:

{Plot3D[Re[qfun[t, K0]], {t, 0, 10}, {K0, 0, 2}, AxesLabel -> {"t", "K0", "Re(Q)"}], 
 Plot3D[Im[qfun[t, K0]], {t, 0, 10}, {K0, 0, 2}, AxesLabel -> {"t", "K0", "Im(Q)"}]}
 // GraphicsRow

Answer 2

I realise this isn't really what you're asking for, but a numerical answer is as close as I've been able to get so far. I'm not convinced there is a closed-form solution for all the non-integer initial conditions Exp[K0 I Pi] (which might explain DSolve's failure). But I'm not convinced there isn't, either.

You could get a numerical solution by building K0 into Q explicitly. That is, make Q = Q[t, K0], and using NDSolveValue (because, as you've discovered, DSolve doesn't seem to be able to do non-integer values of K0).

qfun = NDSolveValue[{D[Q[t, K0], {t, 2}] + Sinh[Q[t, K0]] == 0, 
  (D[Q[t, K0], t] /. t -> 0) == 0, Q[0, K0] == Exp[K0 I Pi]}, 
  Q, {t, 0, 10}, {K0, 0, 2}]

which gives you an InterpolatingFunction in t and K0.

Plotting the real part:

Plot3D[Re[qfun[t, K0]], {t, 0, 10}, {K0, 0, 2}, AxesLabel -> {"t", "K0", "Q"}]

Answer 3

Another option of numerical solution is to make a module and call it for each t and k

f[t0_?NumericQ,k_?NumericQ]:=Module[{sol,t,Q},
   sol=Q/.First@NDSolve[{Q''[t]+Sinh[Q[t]]==0,
                    Q'[0]==0,Q[0]==Exp[k I Pi]},Q,{t,t0,t0+.1}];
   sol[t0]
];

Plot3D[Re[Evaluate[f[t, k]]], {t, 0, 10}, {k, 0, 2}, 
                AxesLabel -> {"t", "K0", "Q"}]

I tried series solution by Maple, but for some reason, the solution it gives does not match the above, no time to find out why. But here it is any way:

restart;
Order:=6:
sol:=dsolve({diff(q(t),t$2) + sinh(q(t)) = 0, 
       D(q)(0)=0, q(0) = exp(k*I*Pi)}, q(t),'series'):

sol:=convert(sol,polynom):
f:=unapply(rhs(sol),t,k);

plot3d(Re(f(t,k)),t=0..10,k=0..2);