# Force NDSolve to choose given branch at a critical point

I am trying to find the solution of f[x,y]==0 by integrating along the curve starting from a point (vars1):

f[x_, y_] = y*(y^2 - x + 1);
vars = {x, y};
varsOt = Through[vars[t]];
vars1 = FindRoot[f[x, 0.5], {x, 0.5}]~Join~{y -> 0.5};
sysDAE0 = {(D[f[x, y], {{x, y}}] /. Thread[vars -> varsOt]).D[varsOt, t] == 0,
Norm[D[varsOt, t]] == 1, (varsOt /. t -> 0) == (vars /. vars1)};
{solx, soly} = NDSolveValue[sysDAE0, vars, {t, -3, 2},
Method -> {"EquationSimplification" -> "Residual"}];
Show[ContourPlot[f[x, y] == 0, {x, -.1, 2}, {y, -2, 2}, PlotPoints -> 100],
ParametricPlot[{solx[t], soly[t]}, {t, -3, 2}, AspectRatio -> 1/2,
PlotStyle -> Red]]


All the root curves are in blue, while the result of the above integration is in red:

Actually, I'd like to choose another branch; instead of take the right part of the $x$-axis, I'd like to catch the left part. To do so, I started by defining a function detecting the critical point (the norm of the gradient of f) in order to trigger WhenEvent. Then, I ask to choose a change the sign of x'[t] when the critical point is detected, hoping it would push the integration to the left of $(1,0)$. But it fails: the integration remains unchanged and still goes to the right of $(1,0)$. Any idea?

critical[x_, y_] = Norm[D[f[x, y], {{x, y}, 1}]];
{solx, soly} =
NDSolveValue[
sysDAE0~Join~{WhenEvent[critical[x[t], y[t]] < 10^-5,
x'[t] -> -x'[t]]}, vars, {t, -3, 2},
Method -> {"EquationSimplification" -> "Residual"}]


Edit To answer some of the comments, here is an example function, which requires the specified Method:

f[x_,y_] = y(9.77516 + 56.827 y^2 + 142.095 y^4 + x^2
(-31.1394 - 162.744 y^2 - 299.61 y^4) + x (9.88476 + 65.3474 y^2 + 180.768 y^4))


Running the NDSolveValue from above returns:

NDSolveValue::ntdv: Cannot solve to find an explicit formula for the derivatives.
Consider using the option Method->{"EquationSimplification"->"Residual"}.


Note that in that case, NDSolve computes the left branch... Then how can I get the right one (in a robust manner)?

Side notes on the equations I chose (answering to xzczd and Chris K's comments): I'd like to solve $f(x,y)=0$. Let's rewrite it $f(X)=0$ with $X=[x\ y]^\top$. Since I wanted to use the efficient NDSolve and the WhenEvent method and get an InterpolatingFunction, I instead looked for $t\mapsto X(t)$ such that $f(X(t))=0$. This implies that $$\dfrac{\mathrm{d}}{\mathrm{d}t}f(X(t))=\nabla_X f(X(t))\cdot X'(t) = 0$$ for some initial conditions $X(0)$ such that $f(X(0))=0$. This gives one (differential) equation on $X$, so the system is undertermined; it can be completed with the following condition: $\|X'(t)\|=1$.

• I have to believe the solution you get has to do with where the approximate step lands at the singular point. MaxStepFraction -> 0.00001 gives a different solution. Jan 30, 2017 at 2:54
• @MichaelE2 Indeed, but how to ensure it goes in a given direction? Do you know what's wrong with my WhenEvent? Jan 30, 2017 at 3:32
• This is a typical problem of a phase transition theory. The solution is degenerated. The degeneration can be removed by an introduction of a small "field" acting "in the direction" of one of the solutions: f[x_, y_] =y*(y^2 - x + 1)-h and assign h to a small value. The closer you are to the bifurcation point the smaller h needs to be. The estimate is h<<Abs[x - 1]^(3/2) (see 1. L. D. Landau and E. M. Lifshitz, Statistical Physics., 3 ed. Pergamon Press, Oxford, 1985, Chapter XIV, Section 144) Jan 30, 2017 at 9:28
• Re WhenEvent: I would think that at each step x'[t] is calculated from the ODEs, so reseting x'[t], if it has any effect at all, has one only for one step. It's also possible that at no step is the norm less than 10^-5 (if the steps were too large). Jan 30, 2017 at 11:11
• @ChrisK Yes; see edit with a more realistic function. Feb 2, 2017 at 2:42

To avoid the ntdv warning (actually in v9.0.1 I only got ndsdtc) and turn to the desired branch we need to

1. modify the 2nd equation a little to avoid Abs and Sqrt.

2. set a more distinct triggering condition for the event.

3. modify the event to x[t] -> x[t] + 10^-2.

Notice I've used NDSolve instead of NDSolveValue because NDSolveValue will only return one solution, which isn't the desired one.

f[x_, y_] =
y (9.77516 + 56.827 y^2 + 142.095 y^4 + x^2 (-31.1394 - 162.744 y^2 - 299.61 y^4) +
x (9.88476 + 65.3474 y^2 + 180.768 y^4));
vars = {x, y};
varsOt = Through[vars[t]];
initx = 1/2; inity = 1/2;
vars1 = FindRoot[f[x, inity], {x, initx}]~Join~{y -> inity};
sysDAE0 = {(D[f[x, y], {{x, y}}] /. Thread[vars -> varsOt]).D[varsOt, t] == 0,
PiecewiseExpand[Norm[D[varsOt, t]], Reals]^2 == 1^2,
(varsOt /. t -> 0) == (vars /. vars1)}
contour = ContourPlot[f[x, y] == 0, {x, -.1, 2}, {y, -2, 2}, PlotPoints -> 100];
critical[x_, y_] = PiecewiseExpand[Norm[D[f[x, y], {{x, y}, 1}]], Reals];
{solx, soly} =
With[{mid = D[critical[x[t], y[t]], t]}, {x, y} /.
NDSolve[{sysDAE0, WhenEvent[mid < 0, x[t] -> x[t] + 10^-2]},
vars, {t, -3, 2}][[2]]];
Show[contour,
ParametricPlot[{solx[t], soly[t]}, {t, -3, 2}, AspectRatio -> 1/2, PlotStyle -> Red]]


• On MMA 11.0, NDSolve returns a function defined over {0.,0.}. Without the WhenEvent, it selects the left part. Now, if I add Print[...] in WhenEvent, as such: WhenEvent[mid < 0, Print["blahblah"]; x[t] -> x[t] + 10^-2]... it selects the right part... (and prints about a hundred blahblah). Feb 2, 2017 at 14:48
• @anderstood Glad to see another bug of v11. (……) Just tested on Cloud, this can be circumvented by adding "DetectionMethod" -> "Sign" in WhenEvent. Feb 2, 2017 at 15:01
• I'm testing the robustness of the method. How would you get the bottom (curved) branch? I changed x to y in x[t] -> x[t] - 10^-2, it returns no error but the branch is not a solution curve, even for different values of the constant. Feb 2, 2017 at 15:22
• @anderstood I'm afraid there won't be a easy (as you expected) way to go to that branch, just observe the result of WhenEvent[mid == 0, {x[t] -> 0.844231, y[t] -> -0.5, "RemoveEvent"}] and you'll understand what I mean. BTW, sorry for my poor math, but why can the corresponding differential equation of a curve be obtained in this way i.e. …(D[f[x, y], {{x, y}}] /. Thread[vars -> varsOt]).D[varsOt, t] == 0, Norm[D[varsOt, t]] == 1…? Feb 3, 2017 at 2:55
• Yes, funny (or hopeless, depending on the point of view :)). Regarding the maths, see bottom of the OP. Feb 3, 2017 at 15:45

Here's a way to numerically get the "curvy forks" of the pitchfork in the second example, using a pseudo-arclength continuation function I hacked together earlier. First, define TrackRootPAL from the edit in this answer and define your f[x,y]. Then get an initial point on the track of interest using FindRoot:

FindRoot[f[0.9, y] == 0, {y, 1.5}]
(* {y -> 0.693003} *)


Then follow it with TrackRootPAL:

tr = TrackRootPAL[{f[x, y]}, {y}, {x, 0.5, 1.1}, 0.9, {0.693003},
NDSolveOpts -> {AccuracyGoal -> 7}]
(* {{y -> InterpolatingFunction[{{0.741048, 1.05351}}, <>]},
{y -> InterpolatingFunction[{{0.741068, 1.05351}}, <>]}} *)


That AccuracyGoal setting took some fiddling to get right.

Finally Plot the results:

Plot[Evaluate[y[x] /. tr], {x, 0, 2}, PlotRange -> {{0, 2}, {-2, 2}}]


You can get the other branch (in this case, trivially y=0) by using a different starting point.

• I think this ends up similar to your approach, but with a different set of equations to solve: sysDAE0 = {f[x[t], y[t]] == 0, Norm[D[varsOt, t]] == 1, (varsOt /. t -> 0) == (vars /. vars1)} instead of sysDAE0 = {(D[f[x, y], {{x, y}}] /. Thread[vars -> varsOt]).D[varsOt, t] == 0, Norm[D[varsOt, t]] == 1, (varsOt /. t -> 0) == (vars /. vars1)}. Why did you set D[f[x, y], {{x, y}}] equal to zero? I'm curious about, but not particularly familiar with this approach, so would appreciate any insights. Feb 2, 2017 at 16:24
• @anderstood oops, forgot to include your handle on my previous comment Feb 2, 2017 at 21:10
• Sorry for my late answer; for the math question, see edit. Regarding your answer: how can you be sure that it will return the curve, and not "bifurcate"? Feb 3, 2017 at 15:43