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]
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
A[0] == 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?"?
Calc
, I think the 5th line of code in the first code block,w0 = w00 - t
, was meant to bew0 = w00 - x
as originally written; I think what's missing is an example value forx
. But thatNDSolve
in the first code block is superfluous. I will update the edit... $\endgroup$