# Plotting Phase Portrait of Duffing Equation

I study a paper that describes a stationary problem where the function $$\Phi_d$$ satisfies the boundary conditions and is governed by a modified Laplace equation with a non-linear term. Here's a breakdown of the equations and the conditions they represent:

1. Equation Details:

• The first equation: $$-k\Delta_d\Phi_d = \omega_d^2(\Phi_n - \beta \Phi_n^3$$

This is a partial differential equation (PDE) where $$\Delta_d$$ denotes the Laplacian in $$d$$ dimensions. The equation includes a linear term proportional to $$\Phi_n$$ and a non-linear term proportional to $$\Phi_n^3$$. The parameter $$k$$ seems to modulate the influence of the spatial derivatives, while $$\omega_d$$ and $$\beta$$ scale the linear and non-linear terms, respectively.

• Boundary conditions: $$\Phi_0 = \Phi_{K+1} = 0$$ These conditions specify that $$\Phi$$ is zero at the boundaries of the domain.

2. Special Case (Non-coupled Particles, $$k=0$$:

• In the limit $$k = 0$$, the spatial derivative terms disappear, simplifying the equation to: $$\ddot{U}_n + \delta\dot{U}_n = \omega_d^2 U_n - \beta \omega_d^2 U_n^3$$. Here, $$U_n$$ seems to represent the state of a system described by the Duffing equation, which is a well-known model for non-linear oscillators with a cubic nonlinearity. The term $$\delta \dot{U}_n$$ introduces damping into the system, representing energy loss.

Here's the Mathematica code for simulating and visualizing the dynamics of the Duffing equation:

(* Parameters of the system *)
ωdSquared = 1;    (* ω_d^2 *)
β = 1;            (* β *)
δ = 0.01;         (* δ, damping coefficient *)

(* Duffing equation setup for use with NDSolve for simulation *)
duffingEquations = {U''[t] + δ U'[t] == ωdSquared U[t] - β U[t]^3,
U[0] == 0.1, U'[0] == 0};

(* Numerical solution of the equations for the trajectory *)
solution = NDSolve[duffingEquations, U, {t, 0, 100}];

(* StreamPlot for phase portrait *)
StreamPlot[{v, ωdSquared u - β u^3 - δ v}, {u, -1.5, 1.5}, {v, -1.5, 1.5}, StreamPoints -> Fine,
StreamColorFunction -> "Rainbow"]

(* Use ParametricPlot for the visualization of the trajectory from the solution *)
ParametricPlot[Evaluate[{U[t], U'[t]} /. solution], {t, 0, 100},
PlotRange -> All, AxesLabel -> {"U(t)", "U'(t)"},
PlotStyle -> Thick]


And of course, it works as I want. My question is the following. Can I create on Mathematica a phase portrait like the one in the paper?

• Could you give a link to the paper? Commented May 5 at 12:47
• @AlexTrounev yes of course! It is in greek thought didaktorika.gr/eadd/handle/10442/45268?locale=en Commented May 5 at 14:26
• Thank you very much. I have AI to translate Dynamics of nonlinear lattice systems: asymptotic behavior and study of the existence and stability of localized oscillations :) Commented May 5 at 14:56
• @AlexTrounev what AI do you use? Commented May 5 at 15:03
• We have several of them on docsie.io :) Commented May 5 at 16:08

The best what we can do with StreamPlot

(*Parameters of the system*)\[Omega]dSquared = 1;    (*\[Omega]_d^2*)
\[Beta] = 1;            (*\[Beta]*)
\[Delta] = 0.5;         (*\[Delta],damping coefficient*)

(*Duffing equation setup for use with NDSolve for simulation*)
duffingEquations = {U''[
t] + \[Delta]  U'[t] == \[Omega]dSquared  U[
t] - \[Beta]  U[t]^3, U[0] == x0, U'[0] == y0};

(*Numerical solution of the equations for the trajectory*)
solution =
ParametricNDSolveValue[
duffingEquations, {U[t], U'[t]}, {t, 0, 100}, {x0, y0}];

Show[StreamPlot[{v, \[Omega]dSquared  u - \[Beta]  u^3 - \[Delta]  \
v}, {u, -2, 2}, {v, -2, 2}, StreamColorFunction -> None,
StreamPoints -> {{{{-1.5, -1.5}, Blue}, {{1.5, 1.5},
Blue}, {{.0255, .0251}, Red}, {{-.05, -.05}, Red}}},
GridLines -> Automatic, PerformanceGoal -> "Quality",
StreamScale -> .15],
ParametricPlot[{Evaluate[ solution[.01, .01]],
Evaluate[ solution[-.01, -.01]], Evaluate[solution[-1.5, -1.5]],
Evaluate[solution[1.5, 1.5]]}, {t, 0, 30}, PlotRange -> All,
AxesLabel -> {"U(t)", "U'(t)"}, PlotStyle -> {Red, Red, Blue, Blue}]]


• (+1) Very good to know that δ = 0.5; Commented May 5 at 14:05
• @cvgmt Thank you. I actually chose this parameter by chance :) Commented May 5 at 15:07
• Using ParametricNDSolve with {a,b} be the initial point.
• We solve {U,U'} simultaneously( see the sol).
• According to @Alex Trounev, should be δ = 0.5;.
Clear["Global*"];
ωdSquared = 1;
β = 1;
δ = 0.5;
duffingEquations = {U''[
t] + δ U'[t] == ωdSquared U[t] - β U[t]^3,
U[0] == a, U'[0] == b};
sol = ParametricNDSolve[
duffingEquations, {U, U'}, {t, 0, 100}, {a, b}];
draw[p_, L_, color_] :=
ParametricPlot[{U[p[[1]], p[[2]]][t], U'[p[[1]], p[[2]]][t]} /.
sol // Evaluate, {t, 0, L}, PlotStyle -> color,
AspectRatio -> Automatic, PerformanceGoal -> "Quality",
PlotRange -> {{-2, 2}, {-2, 2}},
GridLines -> {Range[-2, 2, .5], Range[-2, 2, .5]},
GridLinesStyle -> Directive@{AbsoluteThickness[1], Dotted}];
Show[draw[{-1.5, -1.5}, 25, Directive@{Dashed, Blue}],
draw[{0.01, 0}, 50, Red], draw[{-0.01, 0}, 50, Red],
draw[{1.5, 1.5}, 50, Green]] /.
Line[pts_] :> {AbsoluteThickness[1],
`