Mathematica 9 is now available
Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2008

[Date Index] [Thread Index] [Author Index]

Search the Archive

Complex NDSolve with Multifunctions

  • To: mathgroup at smc.vnet.net
  • Subject: [mg86478] Complex NDSolve with Multifunctions
  • From: "David Park" <djmpark at comcast.net>
  • Date: Wed, 12 Mar 2008 00:13:01 -0500 (EST)

On 4 Mar 2008 Alex Cloninger posted the problem of numerically solving the 
following set of differential equations:

{x'[t]==2 p[t], p'[t]==5/2 I (I x[t])^(3/2), x[0] == 0, p[0] == 1}

He was having difficulty because p'[t] involves a square root, which has two 
values. Mathematica does not return a continuous value of p'[t]. The result 
is a truncated and incomplete solution.

The Presentations package at my web site has a Multivalues capability for 
complex functions. One can have a function that remembers its last value and 
always returns a continuous value by selecting a correct order of roots. 
Using this functionality it is possible to obtain what I believe is a full 
and satisfactory solution to the differential equation and obtain a rather 
detailed picture of what is happening.

Needs["Presentations`Master`"]

The following is a display of a regular Mathematica NDSolve solution. In 
this Manipulate presentation I use two t manipulators, one with a course 
scale and one with a fine scale at the turning point. It is very nice that 
the Manipulate designers allowed that! Of course, when the coarse control is 
used, the fine control is usually out of bounds, but that doesn't matter. 
The red curve is x[t] and the black curve is p[t]. The plot on the right is 
a closeup at one of the turning points. It actually goes the wrong way and 
is why we obtain the truncated solution. The plot on the left is largely 
empty because Mathematica misses the larger solution.

Module[
 {t1 = 0, t2 = 10,
  t3 = 2.75, t4 = 2.85,
  x0 = 0,
  p0 = 1,
  del = .001,
  xpsols},
 Clear[x, p, f, q];
 f[q_?NumberQ] := -(5/2) (-1)^(1/4) q^(3/2);
 xpsols = First@NDSolve[
    {
\!\(\*SuperscriptBox["x", "\[Prime]",
MultilineFunction->None]\)[t] == 2 p[t],

\!\(\*SuperscriptBox["p", "\[Prime]",
MultilineFunction->None]\)[t] == f[x[t]],
     x[0] == x0,
     p[0] == p0}, {x, p}, {t, t1, t2},
    Method -> Automatic,
    MaxSteps -> \[Infinity]];
 {x[t_], p[t_]} = {x[t], p[t]} /. xpsols;
 Manipulate[
  Row[
   {Draw2D[
     {ComplexCurve[p[s], {s, t1, t2},
       PerformanceGoal -> "Quality"],
      Red,
      ComplexCurve[x[s], {s, t1, t2},
       PerformanceGoal -> "Quality"],
      ComplexCirclePoint[p[t], 2, Black, Gray],
      ComplexCirclePoint[x[t], 2, Black, Pink]},
     PlotRange -> {{-100, 175}, {-150, 150}},
     Frame -> True,
     ImageSize -> 250],
    Draw2D[
     {ComplexCurve[p[s], {s, t1, t2},
       PerformanceGoal -> "Quality"],
      Red,
      ComplexCurve[x[s], {s, t1, t2},
       PerformanceGoal -> "Quality"],
      ComplexCirclePoint[p[t], 2, Black, Gray],
      ComplexCirclePoint[x[t], 2, Black, Pink]},
     PlotRange -> 4,
     Frame -> True,
     ImageSize -> 250],
    Draw2D[
     {ComplexCurve[p[s], {s, t3, t4},
       PerformanceGoal -> "Quality"],
      ComplexCirclePoint[p[t], 2, Black, Gray]},
     PlotRange -> {{1 - del, 1 + del}, {.00 - del, .00 + del}},
     Frame -> True,
     ImageSize -> 250]}](* Row *),
  Style["Regular Solution of Differential Equation at Different \
Scales", 16],
  {t, t1, t2, .01, Appearance -> "Labeled"},
  "Fine time control for turning point at t = 2.793",
  {t, t3, t4, .001, Appearance -> "Labeled"}](* Manipulate *)
 ]

The following presents the solution using the Presentations Multivalues 
capability. The full curves are now obtained. Interesting aspects of the 
solution span several different Range scales. The plot on the left shows the 
large scale picture of the trajectories (this is completely missing with the 
regular solution). The middle plot shows the behavior around zero. The right 
hand plot shows the closeup behavior of p[t] at the 'turning point'. It 
looks as if the p[t] curve has a discontinuous slope. If one looks at the 
middle plot one might have expected the solution to go the other way.

Module[
 {t1 = 0, t2 = 10,
  t3 = 2.75, t4 = 2.85,
  x0 = 0,
  p0 = 1,
  del = .001,
  xpsols},
 Clear[x, p, f, q];
 pfunc = Multivalues[
   Null, {-(5/2) (-1)^(1/4) q^(3/2), 5/2 (-1)^(1/4) q^(3/2)}, q];
 f[q_?NumberQ] := Part[CalculateMultivalues[pfunc][q], 1, 1];
 xpsols = First@NDSolve[
    {
\!\(\*SuperscriptBox["x", "\[Prime]",
MultilineFunction->None]\)[t] == 2 p[t],

\!\(\*SuperscriptBox["p", "\[Prime]",
MultilineFunction->None]\)[t] == f[x[t]],
     x[0] == x0,
     p[0] == p0}, {x, p}, {t, t1, t2},
    Method -> Automatic,
    MaxSteps -> \[Infinity]];
 {x[t_], p[t_]} = {x[t], p[t]} /. xpsols;
 Manipulate[
  Row[
   {Draw2D[
     {ComplexCurve[p[s], {s, t1, t2},
       PerformanceGoal -> "Quality"],
      Red,
      ComplexCurve[x[s], {s, t1, t2},
       PerformanceGoal -> "Quality"],
      ComplexCirclePoint[p[t], 2, Black, Gray],
      ComplexCirclePoint[x[t], 2, Black, Pink]},
     PlotRange -> {{-100, 175}, {-150, 150}},
     Frame -> True,
     ImageSize -> 250],
    Draw2D[
     {ComplexCurve[p[s], {s, t1, t2},
       PerformanceGoal -> "Quality"],
      Red,
      ComplexCurve[x[s], {s, t1, t2},
       PerformanceGoal -> "Quality"],
      ComplexCirclePoint[p[t], 2, Black, Gray],
      ComplexCirclePoint[x[t], 2, Black, Pink]},
     PlotRange -> 4,
     Frame -> True,
     ImageSize -> 250],
    Draw2D[
     {ComplexCurve[p[s], {s, t3, t4},
       PerformanceGoal -> "Quality"],
      ComplexCirclePoint[p[t], 2, Black, Gray]},
     PlotRange -> {{1 - del, 1 + del}, {.00 - del, .00 + del}},
     Frame -> True,
     ImageSize -> 250]}](* Row *),
  Style["Differential Equation Solutions at Different Scales", 16],
  {t, t1, t2, .01, Appearance -> "Labeled"},
  "Fine time control for turning point at t = 2.793",
  {t, t3, t4, .001, Appearance -> "Labeled"}](* Manipulate *)
 ]

But that is the limiting case. The last presentation shows the case when 
p[0] is slightly off the critical point at p = 1. Here we see that p[t] 
actually traverses a small loop, the one Alex talked about, and p'[t] is 
quite continuous.

Module[
 {t1 = 0, t2 = 10,
  t3 = 2.790, t4 = 2.81,
  x0 = 0,
  p0 = 1 + .0025 I,
  del = .00002,
  xpsols},
 Clear[x, p, f, q];
 pfunc = Multivalues[
   Null, {-(5/2) (-1)^(1/4) q^(3/2), 5/2 (-1)^(1/4) q^(3/2)}, q];
 f[q_?NumberQ] := Part[CalculateMultivalues[pfunc][q], 1, 1];
 xpsols = First@NDSolve[
    {
\!\(\*SuperscriptBox["x", "\[Prime]",
MultilineFunction->None]\)[t] == 2 p[t],

\!\(\*SuperscriptBox["p", "\[Prime]",
MultilineFunction->None]\)[t] == f[x[t]],
     x[0] == x0,
     p[0] == p0}, {x, p}, {t, t1, t2},
    Method -> Automatic,
    MaxSteps -> \[Infinity]];
 {x[t_], p[t_]} = {x[t], p[t]} /. xpsols;
 Manipulate[
  Row[
   {Draw2D[
     {ComplexCurve[p[s], {s, t1, t2},
       PerformanceGoal -> "Quality"],
      Red,
      ComplexCurve[x[s], {s, t1, t2},
       PerformanceGoal -> "Quality"],
      ComplexCirclePoint[p[t], 2, Black, Gray],
      ComplexCirclePoint[x[t], 2, Black, Pink]},
     PlotRange -> {{-100, 175}, {-150, 150}},
     Frame -> True,
     ImageSize -> 250],

    Draw2D[
     {ComplexCurve[p[s], {s, t1, t2},
       PerformanceGoal -> "Quality"],
      Red,
      ComplexCurve[x[s], {s, t1, t2},
       PerformanceGoal -> "Quality"],
      ComplexCirclePoint[p[t], 2, Black, Gray],
      ComplexCirclePoint[x[t], 2, Black, Pink]},
     PlotRange -> 2,
     Frame -> True,
     ImageSize -> 250],

    Draw2D[
     {ComplexCurve[p[s], {s, t3, t4},
       PerformanceGoal -> "Quality"],
      ComplexCirclePoint[p[t], 2, Black, Gray]},
     PlotRange -> {{1 - del, 1 + del}, {.0025 - del, .0025 + del}},
     Frame -> True,
     ImageSize -> 400]}](* Row *),
  Style["Differential Equation Solutions at Different Scales", 16],
  {t, t1, t2, .01, Appearance -> "Labeled"},
  "Fine time control for turning point at t = 2.799",
  {t, t3, t4, .001, Appearance -> "Labeled"}](* Manipulate *)
 ]



-- 
David Park
djmpark at comcast.net
http://home.comcast.net/~djmpark/




  • Prev by Date: Re: symbolize in v.6
  • Next by Date: FullSimplify interpretation
  • Previous by thread: Re: Adapting function to input type
  • Next by thread: FullSimplify interpretation