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/