1
$\begingroup$

I consider the following system $$ \begin{cases} x'=-x^2 + 3 x - 2 x y=:f(x,y)\\ y'=y^2 - 2 y + x y:=g(x,y) \end{cases} $$

It has a fixed point at $(1,1)$ and the linearization there has a centre, so, in general, the nonlinear system may have either centre or spiral. The trajectories 'contained' in the triangle formed by other fixed points: $(0,0)$, $(3,0)$ and $(0,2)$ (it can be shown rigorously). Numerics also show that indeed this is a nonlinear centre:

f[x_, y_] := -x^2 + 3 x - 2 x y;
g[x_, y_] := y^2 - 2 y + x y;
{xsol, ysol} = NDSolveValue[{x'[t] == f[x[t], y[t]], y'[t] == g[x[t], y[t]], x[0] == 0.5, y[0] == 0.5}, {x, y}, {t, 0, 100}];
{xsol1, ysol1} = NDSolveValue[{x'[t] == f[x[t], y[t]], y'[t] == g[x[t],y[t]], x[0] == 0.8, y[0] == 0.8}, {x, y}, {t, 0, 100}];
ParametricPlot[{{xsol[t], ysol[t]}, {xsol1[t], ysol1[t]}}, {t, 0, 100}, PlotRange -> {0, 2.2}]

enter image description here

I'm looking for an explicit form of these cycled curves. A standard way is to find it at the form $H(x,y)=C=\mathrm{const}$. Then $H(x(t),y(t))=\mathrm{const}$ implies

$$ \frac{\partial H(x,y)}{\partial x} f(x,y)+\frac{\partial H(x,y)}{\partial y}g(x,y)=0. $$

The question is how 'to ask' Mathematica (try to) find such function $H$ in certain form? (It seems that the latter equation can't be handle using DSolve...)

$\endgroup$
4
$\begingroup$

Here's an approach that works for me. There is a timing issue, so I do not claim it is a robust approach. We can use withTimedIntegrate to stop DSolve at the crucial moment. For me, the critical time was around 0.08 sec., so I limited the time allowed for Integrate to 0.01 sec. It yielded the implicit solution the OP seems to be seeking.

ClearAll[withTimedIntegrate];
SetAttributes[withTimedIntegrate, HoldFirst];
withTimedIntegrate[code_, tc_] := 
  Module[{$in}, 
   Internal`InheritedBlock[{Integrate}, Unprotect[Integrate];
    i : Integrate[___] /; ! TrueQ[$in] := 
     Block[{$in = True}, TimeConstrained[i, tc, foo = Inactivate[i, Integrate]]];
    Protect[Integrate];
    code]];

sol1 = withTimedIntegrate[
  DSolve[y'[x] == g[x, y[x]]/f[x, y[x]], y, x], 0.01]

Mathematica graphics

Then we massage the integrals, which luckily evaluate, and adjust some constants (by inspection), and eventually solve for the constant of integration to get $H$:

Simplify[sol1, x > 1] /. HoldPattern@Solve[eq_, _] :> eq /. 
   Inactive[Integrate][e_, u_] :> 
    Inactive[Integrate][e*Dt[u] /. Dt[x] -> 1, x] // Activate;
% /. C[1] -> C[1] (-70)^(2/3)/108 /. y[x] -> y;
H = C[1] - Log[81] /. First@Solve[%, C[1]] // Simplify
(*  2 Log[x] + Log[6 - 2 x - 3 y] + 3 Log[y]  *)

Some of the curves $H = C$:

ContourPlot[H, {x, 0, 3}, {y, 0, 2},
 Contours -> -{0.25, 1.25, 2.5, 4., 5.75, 7.75}, 
 RegionFunction -> 
  Function[{x, y}, x > 0 && y > 0 && 6 - 2 x - 3 y > 0], 
 AspectRatio -> Automatic]

enter image description here

If I limit Integrate to 0.2 sec., I get four solutions. But letting DSolve run without a time-constrained Integrate doesn't yield a result in a time I'm willing to wait.

Of course, $H$ is defined only up to a factor. The integrating factor that leads to the particular $H$ above can be found with the following:

Dt[H]/(g[x, y] Dt[x] - f[x, y] Dt[y]) // Factor
(*  6/(x y (-6 + 2 x + 3 y))  *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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