7
$\begingroup$

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[0] == 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] == 0.7 + 0.5 I
   }, A0[t], {t, 0, tmax}, 
   MaxSteps -> Infinity, AccuracyGoal -> 50, MaxStepSize -> 0.01
  ];

  {Abs[A0[t] /. s2 /. t -> tmax][[1]]}
]
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]

Comparing NDSolve and exact solutions

Solid line – NDSolve solution, dashed line - exact solution

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

A[0] == 0.7 + 0.5*I

enter image description here

A[0] == 10^-6

enter image description here

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?"?

$\endgroup$
6
  • 2
    $\begingroup$ 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. $\endgroup$ Sep 29, 2015 at 21:27
  • $\begingroup$ I just evaluated the code and I ended up with lots of warnings, maybe it needs some cleaning up. $\endgroup$ Sep 29, 2015 at 22:04
  • 1
    $\begingroup$ @RaymondGhaffarianShirazi thanks. I fixed the errors so that the OP's observations are reproducible now. $\endgroup$ Sep 30, 2015 at 0:24
  • $\begingroup$ Thanks, I give it look, I suspect I faced a similar challenge and those sloped vertical lines where produced by plot command actually. $\endgroup$ Sep 30, 2015 at 0:48
  • $\begingroup$ @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... $\endgroup$
    – Michael E2
    Sep 30, 2015 at 10:15

1 Answer 1

10
+100
$\begingroup$

(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[7575] /. {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[7575] /. {u'[t] -> 0, v'[t] -> 0}, {{u[t], v[t]}}] /. equil
(*  {-106153., 192263., 237015.}  *)

StreamPlot[
 ComplexExpand@
  Through[{Re, Im}[
    A0'[t] /. First@Solve[eq[7575], 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]

Mathematica graphics
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] == (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]][[1]]}];

movie = Table[Calc[z];
   ParametricPlot[Through[{Re, Im}[A0[t] /. s2]], {t, 0, 4}, 
    Epilog -> {PointSize[Medium], Red, 
      Point[Through[{Re, Im}[A0[4] /. 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}];

enter image description here

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:

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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