NDSolve's output ignores multiple valid solutions

I'm looking for solutions to a boundary problem involving a non-linear Hamiltonian $$H(q,p) = \frac{1}{4}\left(q^{2}+p^{2}\right)^{2},$$ whose solutions are oscillatory but have a complex time dependence. I'm interested in all possible solutions $\left(q(t),p(t)\right)$ that satisfy the following boundary conditions:

$$\begin{cases} q(0)&=-1 \\ q(\pi)&=1 \end{cases}$$

and I am absolutely sure there are a lot (maybe an infinity) of them. When I ask Mathematica to solve the boundary problem

NDSolve[{q'[t] == p[t] (p[t]^2 + q[t]^2),
p'[t] == -q[t] (p[t]^2 + q[t]^2), q == -1, q[Pi] == 1}, q[t], {t, 0, Pi}]

I get only one solution, which looks like and satisfies the boundary problem. What I can't figure out is that, manually, I found another solution:

NDSolve[{q'[t] == p[t] (p[t]^2 + q[t]^2),
p'[t] == -q[t] (p[t]^2 + q[t]^2), q == -1, p == 1.200859},
q[t], {t, 0, Pi}]

whose graph is

How can I manipulate NDSolve such that it displays more solutions? Since they may be infinite, not all can be displayed, but why is Mathematica just choosing a particular solution in a set of infinite ones?

Update

You seem correct QuantumBrick that the Shooting method is better:

sols = Map[First[
NDSolve[{q'[t] == p[t] (p[t]^2 + q[t]^2),
p'[t] == -q[t] (p[t]^2 + q[t]^2),
q == -1, q[Pi] == 1}, {q, p}, {t, 0, Pi},
Method ->  "BoundaryValues" -> {"Shooting",
"StartingInitialConditions" -> {p == #}}]] &, Range[0.25, 2, 0.25]];

Plot[Evaluate[q[t] /. sols], {t, 0, Pi}] Introducing small error into the starting conditions to find other approximate answers (which is similar to your manual answer)

sol = Table[
NDSolve[{q'[t] == p[t] (p[t]^2 + q[t]^2),
p'[t] == -q[t] (p[t]^2 + q[t]^2),
q == -RandomReal[{0.99, 1.01}],
q[Pi] == RandomReal[{0.99, 1.01}]}, q, {t, 0, Pi}], {10}];

Plot[Table[q[t] /. sol[[i]], {i, 1, 10}], {t, 0, Pi}] • Thansk for the answer, but I sort of solved it right this moment using the Shooting method with starting conditions on the momenta. Is your method more precise then this one? – QuantumBrick Aug 1 '16 at 16:49
• @QuantumBrick The Shooting method may be better ... – Young Aug 1 '16 at 16:54
• "I sort of solved it right this moment" - next time, @Quantum Brick, consider posting an answer to your own question if you manage to figure out a solution after asking it. :) – J. M. will be back soon Aug 1 '16 at 17:08
• @J.M. Should I just delete my answer? – Young Aug 1 '16 at 17:10
• @Young, please don't, as you've written it already. I was just suggesting to the OP that answering your own question is kosher here, should s/he be in that situation the next time. – J. M. will be back soon Aug 1 '16 at 17:11

The ode can be solved symbolically, except DSolve runs into trouble with the branches of Sqrt[] and we end up with a general solution that is essentially -Abs[solution]. As result, DSolve[] can't solve the boundary conditions (they cannot be satisfied since the computed q is nonpositive for all initial conditions). But all is not lost.

ode = {q'[t] == p[t] (p[t]^2 + q[t]^2), p'[t] == -q[t] (p[t]^2 + q[t]^2)};
bcs = {q == -1(*,q[Pi] == 1*)}; (* change to half an IVP *)
dsols = DSolve[{ode, bcs}, {p, q}, t];
(* DSolve/Solve warnings about inverse function being used *)

DSolve returns two solutions (with the integration constant C). We can, it turns out if we check, transform one of the solutions into a differentiable solution. Here we define the solution as a pair of functions qfn and pfn and check that they solve the ode.

With[{sol = Last@dsols},              (* the Last leads to a solution *)
{qfn = Evaluate[PowerExpand@Simplify[q[t] /. sol] /. t -> #] &,
pfn = Evaluate[PowerExpand@Simplify[p[t] /. sol] /. t -> #] &}]
With[{sol = First@dsols},             (* the First does not lead to a solution *)
{qfn2 = Evaluate[PowerExpand@Simplify[q[t] /. sol] /. t -> #] &,
pfn2 = Evaluate[PowerExpand@Simplify[p[t] /. sol] /. t -> #] &}];
(*
{-Cos[C + #1 + #1 Tan[C]^2] Sec[C] &,
Sec[C] Sin[C + #1 + #1 Tan[C]^2] &}
*)

ode /. {q -> qfn, p -> pfn} // Simplify
ode /. {q -> qfn2, p -> pfn2} // Simplify  (* does not satisfy ode *)
(*
{True, True}

{Sec[C] Sin[t + C + t Tan[C]^2] == 0,
Cos[t + C + t Tan[C]^2] Sec[C] == 0}
*)

We cannot solve the boundary condition q[Pi] == 1 symbolically, but NSolve[] can handle it, if we restrict the domain. First, let's look at what we're going to solve:

eq = q[Pi] - 1 /. Last@dsols // Simplify // PowerExpand;

Plot[eq /. C -> c2, {c2, 0, 1.4}] NSolve won't be able to solve the equation in a neighborhood of Pi/2, since there are infinitely many solutions. It also has trouble with the solution C -> 0 for some reason.

Here we compute 1001 solutions:

bcsols = Join[
{{C -> 0}},
NSolve[{eq == 0, 0 < C < Pi/2 - 0.01}, C]
];
Length@bcsols
(*  1001  *)

Here's a look at the first twenty:

Plot[Take[qfn[t] /. bcsols, 20] // Evaluate, {t, 0, Pi}] Here's a check of the boundary conditions of the 1001-st solution:

qfn /@ {0, Pi} /. bcsols[]
(*  {-1, 1.}  *)

This will get you another 9,999 solutions on the other side of C == 0:

bcsols2 = NSolve[{eq == 0, -Pi/2 + 0.01 < C < -0.01}, C];
Length@bcsols2
(*  9999  *)

Not NSolve is a bit finicky: The constraint C < 0 is not good enough; you need C less than a (not too small) negative number.

• This is way more than what I've asked. I'll study your answer in a few moments, but thank you very much for the insights you provided. This is an exceptional answer. – QuantumBrick Aug 1 '16 at 19:11
• @QuantumBrick You're welcome...& thanks! – Michael E2 Aug 1 '16 at 19:21
• @QuantumBrick Actually you seemed to hint that you were looking for "a lot (maybe an infinity) of them" -- ha-ha :) – Michael E2 Aug 1 '16 at 19:27
• (For some reason I can't tag you) - Physics told me there should be an infinity of them, but the algebraic approach you performed was quite unexpected and welcome ;) – QuantumBrick Aug 1 '16 at 19:45
• @QuantumBrick Authors of posts are automatically tagged without the @ stuff. From a quick math glance, too, one expects infinitely many. – Michael E2 Aug 1 '16 at 19:48

An alternative, perhaps simpler, approach makes use of the fact that the Hamiltonian is a constant of the motion here. This can be validated by constructing the time-derivative of the Hamiltonian.

eqs = {q'[t] == p[t] (p[t]^2 + q[t]^2), p'[t] == -q[t] (p[t]^2 + q[t]^2)};
q[t] eqs[[1, 2]] + p[t] eqs[[2, 2]]
(* 0 *)

Setting the Hamiltonian equal to ω^2 then yields

eqs1 = eqs /. (p[t]^2 + q[t]^2) -> ω
(* {Derivative[q][t] == ω p[t], Derivative[p][t] == -ω q[t]} *)

which DSolve handles without difficulty.

s = Collect[DSolve[{eqs1, q == -1, q[Pi] == 1}, {p[t], q[t]}, t],
{Sin[t ω], Cos[t ω]}, Simplify] // Flatten
(* {p[t] -> Cos[t ω] Cot[(π ω)/2] + Sin[t ω], q[t] -> -Cos[t ω] + Cot[(π ω)/2] Sin[t ω]} *)

Edit: Finally, ω is determined by p[t]^2 + q[t]^2 == ω

Total[#^2 & /@ ({p[t], q[t]} /. s)] // FullSimplify
(* Csc[(π ω)/2]^2 *)

So, the eigenvalue equation is

ω Sin[(π ω)/2]^2 - 1 == 0

Plotting this function indicates the locations of the roots,

Plot[ω Sin[(π ω)/2]^2 - 1, {ω, 0, 12}] And the roots themselves are given by

freq = ω /. NSolve[ω Sin[(π ω)/2]^2 - 1 == 0 && .4 < ω < 20, ω]
(* {1., 1.33333, 2.44206, 3.64927, 4.31956, 5.72552, 6.26172, 7.76635,
8.22672, 9.79293, 10.2027, 11.812, 12.185, 13.8267, 14.1712, 15.8383,
16.16, 17.8479, 18.1508, 19.8559} *)

p1 = Plot[Evaluate[q[t] /. s /. ω -> freq[[1 ;; 6]]], {t, 0, Pi}] The derivation above misses some solutions. Apply DSolve to eqs1 without the boundary conditions.

{q[t], p[t]} /. DSolve[{eqs1}, {p[t], q[t]}, t] // First
(* {C Cos[t ω] + C Sin[t ω], C Cos[t ω] - C Sin[t ω]} *)

The first boundary condition yields

% /. t -> 0
(* {C, C} *)

Consequently, C == -1 and ω == 1 + C^2. The second boundary condition then yields

Reduce[(%%[] /. {C -> -1, t -> Pi}) == 1, C] // FullSimplify
(* (Sin[π ω] == 0 && Cos[π ω] == -1) || (Sin[π ω] != 0 && C == Cot[(π ω)/2]) *)

The second result is the one obtained earlier. The first, however, is new. It is satisfied by ω any odd integer. Corresponding solutions are, for instance,

p2 = Plot[Evaluate[{-Cos[t ω] + Sqrt[ω - 1] Sin[t ω], -Cos[t ω] - Sqrt[ω - 1] Sin[t ω]} /.
ω -> {3, 5}], {t, 0, Pi}, PlotStyle -> Dashed] Between them, plots p1 and p2 depict all solutions for ω < 6.

Show[p1, p2] Incidentally, one might have expected DSolve with the boundary conditions and

SetOptions[Solve, Method -> Reduce];

to return both sets of solutions but it does not.

• I'm really curious about why DSolve couldn't map the second solution. I think this might be a good question on its own. – QuantumBrick Aug 15 '16 at 17:47
• I also believe some other solutions are missing. You chose $C=\sqrt{\omega-1}$, but $C=-\sqrt{\omega-1}$ should also be a solution. – QuantumBrick Aug 15 '16 at 17:59
• @QuantumBrick p2 shows curves for {-Cos[t ω] + Sqrt[ω - 1] Sin[t ω], -Cos[t ω] - Sqrt[ω - 1] Sin[t ω]}. The first expression in the list evaluates the positive square root, and the second expression the negative square root. – bbgodfrey Aug 15 '16 at 20:27
• @QuantumBrick With respect to your first comment, the second solution is, in essence, the solution of a Sturm-Liouville problem. The related function, DEigensystem cannot treat inhomogeneous boundary conditions, and it seems likely that DSolve cannot either. Certainly, all Sturm-Liouville problems in the DSolve documentation have homogeneous boundary conditions. If you were to ask this as a new question, the likely answer would be that Wolfram, Inc has not yet implemented this capability. – bbgodfrey Aug 15 '16 at 21:13