12
$\begingroup$

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

enter image description here

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

enter image description here .

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?

$\endgroup$
8
$\begingroup$

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}]

enter image description here


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}]

enter image description here

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – QuantumBrick Aug 1 '16 at 16:49
  • 1
    $\begingroup$ @QuantumBrick The Shooting method may be better ... $\endgroup$ – Young Aug 1 '16 at 16:54
  • $\begingroup$ "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. :) $\endgroup$ – J. M. will be back soon Aug 1 '16 at 17:08
  • $\begingroup$ @J.M. Should I just delete my answer? $\endgroup$ – Young Aug 1 '16 at 17:10
  • $\begingroup$ @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. $\endgroup$ – J. M. will be back soon Aug 1 '16 at 17:11
9
$\begingroup$

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}]

Mathematica graphics

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}]

Mathematica graphics

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – QuantumBrick Aug 1 '16 at 19:11
  • $\begingroup$ @QuantumBrick You're welcome...& thanks! $\endgroup$ – Michael E2 Aug 1 '16 at 19:21
  • 1
    $\begingroup$ @QuantumBrick Actually you seemed to hint that you were looking for "a lot (maybe an infinity) of them" -- ha-ha :) $\endgroup$ – Michael E2 Aug 1 '16 at 19:27
  • $\begingroup$ (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 ;) $\endgroup$ – QuantumBrick Aug 1 '16 at 19:45
  • $\begingroup$ @QuantumBrick Authors of posts are automatically tagged without the @ stuff. From a quick math glance, too, one expects infinitely many. $\endgroup$ – Michael E2 Aug 1 '16 at 19:48
6
$\begingroup$

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}]

enter image description here

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}]

enter image description here

Addendum

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]

enter image description here

Between them, plots p1 and p2 depict all solutions for ω < 6.

Show[p1, p2]

enter image description here

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.

$\endgroup$
  • $\begingroup$ I'm really curious about why DSolve couldn't map the second solution. I think this might be a good question on its own. $\endgroup$ – QuantumBrick Aug 15 '16 at 17:47
  • $\begingroup$ 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. $\endgroup$ – QuantumBrick Aug 15 '16 at 17:59
  • $\begingroup$ @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. $\endgroup$ – bbgodfrey Aug 15 '16 at 20:27
  • 1
    $\begingroup$ @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. $\endgroup$ – bbgodfrey Aug 15 '16 at 21:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.