Phase Portraits of Nonlinear Systems

Nonlinear Autonomous Systems

A planar nonlinear autonomous system has the form \(\dot{x}=f(x,y)\), \(\dot{y}=g(x,y)\) where \(f\) and \(g\) are arbitrary (possibly nonlinear) functions. Unlike linear systems, the qualitative behaviour can vary across the phase plane: a single system can have multiple equilibria, limit cycles, homoclinic orbits, and chaotic regions in different parts of the plane.

Equilibria and Linearisation

An equilibrium is a point \((x^*,y^*)\) satisfying \(f(x^*,y^*)=g(x^*,y^*)=0\). Near each equilibrium, the system behaves like a linear system governed by the Jacobian matrix \(J = \partial(f,g)/\partial(x,y)\). The eigenvalues of \(J\) at the equilibrium determine the local classification — saddle, node, spiral, or center. Away from equilibria, the nonlinear terms dominate and the global portrait can differ substantially from any local linearisation.

Enable Equilibria to locate and classify the fixed points automatically: Newton's method is seeded across the domain to solve \(f=g=0\), and each root is classified from the trace \(\tau\) and determinant \(\Delta\) of the numerical Jacobian. Markers follow the convention filled disc = stable (node or spiral), hollow disc = unstable, × = saddle (\(\Delta<0\)), and ring with a dot = centre (\(\tau=0,\ \Delta>0\)); a legend lists the types present. A linear centre is borderline — the true nonlinear behaviour there may be a slow spiral.

Nullclines

The x-nullcline is the set of points where \(\dot{x}=f(x,y)=0\) (the field is vertical there) and the y-nullcline is where \(\dot{y}=g(x,y)=0\) (the field is horizontal). Equilibria are exactly the points where the two nullclines intersect. Enable Nullclines to overlay them — the x-nullcline is drawn solid and the y-nullcline dashed, in your chosen colour — and use the regions they carve up to reason about the sign of \(\dot{x}\) and \(\dot{y}\) and the qualitative flow.

Limit Cycles

A limit cycle is an isolated closed orbit that neighbouring trajectories spiral toward (stable) or away from (unstable). Limit cycles cannot occur in linear systems. The Van der Pol and FitzHugh-Nagumo systems in the example dropdown both exhibit stable limit cycles. The Poincaré-Bendixson theorem guarantees that if a trajectory is bounded and avoids equilibria, it must approach a limit cycle.

Numerical Solvers

• Tsit5 — Tsitouras 5(4): A 7-stage explicit Runge-Kutta pair (Tsitouras, 2011) with First-Same-As-Last (FSAL) property, giving effectively 6 function evaluations per accepted step. Provides 5th-order accuracy with an embedded 4th-order error estimator for adaptive step control. The default solver in Julia's DifferentialEquations.jl for non-stiff problems. Default choice; recommended for most systems.
• Dormand-Prince 5(4): The classic adaptive 5(4) pair (Dormand & Prince, 1980) used in MATLAB's ode45 and SciPy's RK45. Also FSAL with 7 stages. Very reliable for general non-stiff systems.
• DOP853 — Dormand-Prince 8(5,3): A high-order explicit Runge-Kutta method of order 8 with embedded 5th and 3rd-order error estimators (Hairer, Nørsett & Wanner). Uses 12 stages per step but achieves much higher accuracy per step than 5th-order methods. Recommended when high precision is needed (fine/very fine step size settings). Implementation follows the SciPy/Hairer reference.
• RK4 — Classical Runge-Kutta: The standard 4th-order fixed-step Runge-Kutta method. Uses 4 stages with no adaptive error control — the step size is determined entirely by the Step size selector. Good for pedagogical comparison and systems where uniform sampling is desired.
• SDIRK2 (Singly Diagonally Implicit RK, 2nd order): An A-stable, L-stable implicit solver designed for stiff systems. Uses Newton iterations at each stage with numerical Jacobian and adaptive step size via embedded error estimation. Recommended for Van der Pol with large \(\mu\), Brusselator, and Selkov models.
• Rodas5P — Rosenbrock-Wanner 5(4): An 8-stage L-stable Rosenbrock method (Steinebach, 2023) with embedded 4th-order error estimator. Linearly implicit — no Newton iteration; one 2×2 W-matrix factorisation \(W = I/(\gamma h) - J\) is reused across all 8 stages. Typically the fastest solver here for stiff non-linear systems. Coefficients ported from OrdinaryDiffEq.jl.
• RadauIIA5 — Implicit Runge-Kutta 5(3): A 3-stage fully-implicit method (Hairer-Wanner Vol II §IV.8) based on Radau quadrature. A- and L-stable, stiffly accurate. Uses simplified Newton on the 6×6 stage system with the Jacobian frozen at \(u_n\). The reference stiff solver in scientific computing.

Expression Syntax

• Use standard math notation: x^2, sin(x), exp(-x), sqrt(y), pi, e.
• Multiplication must be explicit: x*y not xy.
• Powered by math.js — see its documentation for the full list of supported functions.

How to Use

• Enter \(f(x,y)\) and \(g(x,y)\) in the expression fields, or select an example from the dropdown.
• Set the domain \([x_\min, x_\max]\times[y_\min, y_\max]\) and adjust arrow settings.
• Click Generate to draw the direction field. Then click anywhere on the phase plane to integrate a trajectory through that initial condition, or type an exact \((x_0, y_0)\) in the Initial condition box and press Add trajectory (keyboard-friendly).
• Tsit5 (default) handles most systems. For high-precision work, try DOP853. For stiff problems (e.g. Van der Pol with large \(\mu\)), Rodas5P is usually fastest, followed by RadauIIA5 and SDIRK2. RK4 provides fixed-step integration for comparison. The Step size selector controls the error tolerance (or fixed step for RK4).
• Open the Time Series panel to plot \(x(t)\) and/or \(y(t)\) for all current trajectories.
• Use Download PNG to save the current portrait.