# Jumps in NDSolve results

I need to compute using NDSolve routine, some function $F(x)$, having two possible values $F_1(x)$ and $F_2(x)$ depending on whether the argument exceeds some critical value $x>x_c$. The problem is that NDSolve routine return "jumping results" for the region $x>x_c$.

As an example of a problem assume following system (some simple forced pendulum with nonlinear frequncy shift, so-called "fold over" effect):

w00 = 7000;
G0 = 50;
Q = 14000;
hU = 0.6*G0;
w0 = w00 - x;

(* example NDSolve call with sample value for x -- see below for actual usage *)
Block[{x = 8000},
NDSolve[
{A'[t] + G0*A[t] + I*(w0+Q*Abs[A[t]]^2)*A[t] == I*hU, A == 10^-6},
A[t], {t, 0, tmax},
MaxSteps -> Infinity, AccuracyGoal -> 50, MaxStepSize -> 0.01
];
]


I am interested in $|A(t)|^2$ dynamics depending on the parameter x. Here $G0,Q,hU$ are some constants

Exact solution:

exact = ContourPlot[
A == (hU^2/((G0)^2 + (w0 + Q*A)^2)),
{t, 6000, 9000}, {A, 0, 0.2},
PlotPoints -> 40, ContourStyle -> {Dashed, Thick},
PlotRange -> All
]


Numerical solution:

Calc[x_] := Module[{
w0 = 7000, G0 = 50, Q = 14000, At = 10^-6, tmax = 4,
hU, w00, G1
},

hU = 0.6*G0;
w00 = w0 - x + Q*Abs[A0[t]]^2;
G1 = 1.0 G0;

s2 = NDSolve[{
A0'[t] + G0*A0[t] + I*w00*A0[t] == I*hU,
A0 == 0.7 + 0.5 I
}, A0[t], {t, 0, tmax},
MaxSteps -> Infinity, AccuracyGoal -> 50, MaxStepSize -> 0.01
];

{Abs[A0[t] /. s2 /. t -> tmax][]}
]
Dynamic[z]
For[h1 = {}; z = 9000, z > 6000, z = z - 50,
AppendTo[h1, {z, First@Calc[z]^2}];
]
numerical = ListPlot[h1, PlotRange -> All, Joined -> True, PlotStyle -> Thick]


Comparing both:

Show[exact, numerical] Solid line – NDSolve solution, dashed line - exact solution

The result and number of jumps greatly depends on initial conditions:

A == 0.7 + 0.5*I A == 10^-6 The problem is I want to divide the two "branches" of the function $F_1$ and $F_2$ and plot them independently. How can I do this using NDSolve and force it to go along the chosen branch and not "jump?"?

• To closers: yes, this is an old, unanswered question. But (a) it is upvoted, (b) OP has been seen only 2 months ago, and (c) I cannot agree that it is at all unclear. Potentially it is too localized, but I'm not sure about that. I appreciate the help but this is not quite what I had in mind by closing old questions. I think the best thing to do is to bounty it--please leave a comment if you think it should still be closed. – Oleksandr R. Sep 29 '15 at 21:27
• I just evaluated the code and I ended up with lots of warnings, maybe it needs some cleaning up. – Raymond Ghaffarian Shirazi Sep 29 '15 at 22:04
• @RaymondGhaffarianShirazi thanks. I fixed the errors so that the OP's observations are reproducible now. – Oleksandr R. Sep 30 '15 at 0:24
• Thanks, I give it look, I suspect I faced a similar challenge and those sloped vertical lines where produced by plot command actually. – Raymond Ghaffarian Shirazi Sep 30 '15 at 0:48
• @OleksandrR. Judging from the definition of Calc, I think the 5th line of code in the first code block, w0 = w00 - t, was meant to be w0 = w00 - x as originally written; I think what's missing is an example value for x. But that NDSolve in the first code block is superfluous. I will update the edit... – Michael E2 Sep 30 '15 at 10:15

(Update: I forgot to copy some of the code)

For values of z = x beyond the critical value (just below z == 7500), there are three real equilibria, two stable spirals and one (unstable) saddle.

eq[x_] := Module[{q = x},
w0 = 7000;
G0 = 50;
Q = 14000;
hU = 0.6*G0;
w00 = w0 - q + Q*Abs[A0[t]]^2;
At = 10^-6;
G1 = 1.0 G0;
tmax = 4;
e1 = D[A0[t], t] + G0*A0[t] + I*w00*A0[t] == I*hU];

realeq[x_] := (* the real and imaginary parts in terms of A0[t] = u[t] + I v[t] *)
ComplexExpand@
Through[{Re, Im}[Subtract @@ eq[x] /. {A0 -> (u[#] + I v[#] &)}]] == {0, 0}

equil = NSolve[realeq /. {u'[t] -> 0, v'[t] -> 0}, {u[t], v[t]}, Reals]
(*
{{u[t] -> -0.163309,  v[t] -> 0.0483454},
{u[t] ->  0.207442,  v[t] -> 0.0832793},
{u[t] -> -0.0560372, v[t] -> 0.00528008}}
*)


Stability:

Det@D[First@realeq /. {u'[t] -> 0, v'[t] -> 0}, {{u[t], v[t]}}] /. equil
(*  {-106153., 192263., 237015.}  *)

StreamPlot[
ComplexExpand@
Through[{Re, Im}[
A0'[t] /. First@Solve[eq, A0'[t]] /. {A0[t] -> u + I v}]],
{u, -0.3, 0.3}, {v, -0.3, 0.3},
Epilog -> {PointSize[Medium], Green, Point[{u[t], v[t]} /. equil]},
StreamPoints -> Fine] Fig. 1. Phase plot with separatrices showing the basins of attraction.

As z changes, the basins of attraction for each spiral point change, and as they move, the initial conditions changes which basin it lies in.

Calc[x_] := Module[{q = x},
w0 = 7000;
G0 = 50;
Q = 14000;
hU = 0.6*G0;
w00 = w0 - q + Q*Abs[A0[t]]^2;
At = 10^-6;
G1 = 1.0 G0;
tmax = 4;
e1 = D[A0[t], t] + G0*A0[t] + I*w00*A0[t] == I*hU;
s2 = NDSolve[{e1, A0 == (0.7 + 0.5 I) + 0*10^-6}, {A0}, {t, 0,
tmax}, MaxSteps -> Infinity, AccuracyGoal -> 50,
MaxStepSize -> 0.01];
{Evaluate[Abs[A0[t] /. s2 /. t -> tmax]][]}];

movie = Table[Calc[z];
ParametricPlot[Through[{Re, Im}[A0[t] /. s2]], {t, 0, 4},
Epilog -> {PointSize[Medium], Red,
Point[Through[{Re, Im}[A0 /. First@s2]]], Green,
Point[{u[t], v[t]} /.
NSolve[ComplexExpand@
Through[{Re, Im}[
Subtract @@ e1 /. {A0[t] -> u[t] + I v[t], A0'[t] -> 0}]] == {0, 0},
{u[t], v[t]}, Reals]]},
PlotRange -> 1, MaxRecursion -> 15,
PlotLabel -> HoldForm["z"] == z],
{z, 7550, 7650}]; To get a result like the one the OP wants, one would need to adjust the initial condition as a function of z. But if, as it seems, the OP just wants to get a function representing the curve in the contour plot, one could use a different method, such as the ones I used here: