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/