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$.
MaxStepFraction -> 0.00001
gives a different solution. $\endgroup$ – Michael E2 Jan 30 '17 at 2:54WhenEvent
? $\endgroup$ – anderstood Jan 30 '17 at 3:32f[x_, y_] =y*(y^2 - x + 1)-h
and assignh
to a small value. The closer you are to the bifurcation point the smallerh
needs to be. The estimate ish<<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) $\endgroup$ – Alexei Boulbitch Jan 30 '17 at 9:28WhenEvent
: I would think that at each stepx'[t]
is calculated from the ODEs, so resetingx'[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 than10^-5
(if the steps were too large). $\endgroup$ – Michael E2 Jan 30 '17 at 11:11