# 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[0] == -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[0] == -1, p[0] == 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?

• Have you seen "SymplecticPartitionedRungeKutta" Method for NDSolve? Commented Aug 1, 2016 at 19:56
• @MichaelE2 Man... I've never heard of that, and it's awesome. If you send me one more piece of info I'll be happy to cite you on the paper I'm working at ;) Commented Aug 1, 2016 at 19:58
• Very glad to help. :) Commented Aug 1, 2016 at 20:06
• I was just reminded that there is a bug in SPRK (I ran into it). I posted a fix: mathematica.stackexchange.com/questions/90066/… Commented Aug 7, 2016 at 16:34

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[0] == -1, q[Pi] == 1}, {q, p}, {t, 0, Pi},
Method ->  "BoundaryValues" -> {"Shooting",
"StartingInitialConditions" -> {p[0] == #}}]] &, 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[0] == -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? Commented Aug 1, 2016 at 16:49
• @QuantumBrick The Shooting method may be better ... Commented Aug 1, 2016 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. :) Commented Aug 1, 2016 at 17:08
• @J.M. Should I just delete my answer? Commented Aug 1, 2016 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. Commented Aug 1, 2016 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[0] == -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[2]). 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[2] + #1 + #1 Tan[C[2]]^2] Sec[C[2]] &,
Sec[C[2]] Sin[C[2] + #1 + #1 Tan[C[2]]^2] &}
*)

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

{Sec[C[2]] Sin[t + C[2] + t Tan[C[2]]^2] == 0,
Cos[t + C[2] + t Tan[C[2]]^2] Sec[C[2]] == 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[2] -> 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[2] -> 0 for some reason.

Here we compute 1001 solutions:

bcsols = Join[
{{C[2] -> 0}},
NSolve[{eq == 0, 0 < C[2] < Pi/2 - 0.01}, C[2]]
];
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[[1001]]
(*  {-1, 1.}  *)


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

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


Not NSolve is a bit finicky: The constraint C[2] < 0 is not good enough; you need C[2] 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. Commented Aug 1, 2016 at 19:11
• @QuantumBrick You're welcome...& thanks! Commented Aug 1, 2016 at 19:21
• @QuantumBrick Actually you seemed to hint that you were looking for "a lot (maybe an infinity) of them" -- ha-ha :) Commented Aug 1, 2016 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 ;) Commented Aug 1, 2016 at 19:45
• @QuantumBrick Authors of posts are automatically tagged without the @ stuff. From a quick math glance, too, one expects infinitely many. Commented Aug 1, 2016 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[1][q][t] == ω p[t], Derivative[1][p][t] == -ω q[t]} *)


which DSolve handles without difficulty.

s = Collect[DSolve[{eqs1, q[0] == -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[2] Cos[t ω] + C[1] Sin[t ω], C[1] Cos[t ω] - C[2] Sin[t ω]} *)


The first boundary condition yields

% /. t -> 0
(* {C[2], C[1]} *)


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

Reduce[(%%[[1]] /. {C[2] -> -1, t -> Pi}) == 1, C[1]] // FullSimplify
(* (Sin[π ω] == 0 && Cos[π ω] == -1) || (Sin[π ω] != 0 && C[1] == 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. Commented Aug 15, 2016 at 17:47
• I also believe some other solutions are missing. You chose $C[1]=\sqrt{\omega-1}$, but $C[1]=-\sqrt{\omega-1}$ should also be a solution. Commented Aug 15, 2016 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. Commented Aug 15, 2016 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. Commented Aug 15, 2016 at 21:13