Plotting Phase Portrait of Duffing Equation

I study a paper that describes a stationary problem where the function 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:

This is a partial differential equation (PDE) where denotes the Laplacian in dimensions. The equation includes a linear term proportional to and a non-linear term proportional to . The parameter seems to modulate the influence of the spatial derivatives, while and scale the linear and non-linear terms, respectively.

• Boundary conditions:

These conditions specify that $ \Phi $ is zero at the boundaries of the domain.

2. Special Case (Non-coupled Particles, :)

• In the limit $ k = 0 $, the spatial derivative terms disappear, simplifying the equation to:

Here, 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 introduces damping into the system, representing energy loss.

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

% Parameters
omega_d_squared = 1;
beta = 1;
delta = 0.5;

% Duffing equations
duffingEquations = @(t, y) [y(2); omega_d_squared * y(1) - beta * y(1)^3 - delta * y(2)];

% Initial conditions and time span
initial_conditions = [
    -1.5, -1.5;
    0.01, 0;
    -0.01, 0;
    1.5, 1.5
];
time_span = [0, 50];

% Colors for the plots
colors = [
    0, 0, 1;    % Blue
    1, 0, 0;    % Red
    1, 0, 0;    % Red
    0, 1, 0     % Green
];

% Line styles
line_styles = {'--', '-', '-', '-'};

% Plot the solutions
figure;
hold on;
box on; % Add box around the plot

for i = 1:size(initial_conditions, 1)
    [t, sol] = ode45(duffingEquations, time_span, initial_conditions(i, :));
    plot(sol(:, 1), sol(:, 2), 'Color', colors(i, :), 'LineStyle', line_styles{i}, 'LineWidth', 1);
end

% Customize the plot
axis([-2 2 -2 2]);
grid on;
grid minor;
set(gca, 'XTick', -2:0.5:2, 'YTick', -2:0.5:2, 'GridLineStyle', ':');
set(gca, 'Box', 'on'); 



% Labels and title
xlabel('$U_n$', 'Interpreter', 'latex');
ylabel('$U_n^\prime$', 'Interpreter', 'latex'); 
title('Duffing Oscillator Phase Portrait');
text(1, 1.7, '\omega_d^2 = 1', 'FontSize', 10);
text(1, 1.5, '\beta = 1', 'FontSize', 10);
hold off;

My question is the following, how could I those arrows on the trajectories

[Sam Chak]

Hi, Have you resolved the issue using quiver? Have you tried the annotation command if the quiver's arrowhead too small in the print?

[Athanasios Paraskevopoulos]

Hello@Sam Chak I tried the quiver and gave the following! The arrowheads are really to small for a paper. Now for annotation I have never used before

% Parameters
omega_d_squared = 1;
beta = 1;
delta = 0.5;

% Duffing equations
duffingEquations = @(t, y) [y(2); omega_d_squared * y(1) - beta * y(1)^3 - delta * y(2)];

% Initial conditions and time span
initial_conditions = [
    -1.5, -1.5;
    0.01, 0;
    -0.01, 0;
    1.5, 1.5
];
time_span = [0, 50];

% Colors for the plots
colors = [
    0, 0, 1;    % Blue
    1, 0, 0;    % Red
    1, 0, 0;    % Red
    0, 1, 0     % Green
];

% Line styles
line_styles = {'--', '-', '-', '-'};

% Plot the solutions
figure;
hold on;
box on; % Add box around the plot

for i = 1:size(initial_conditions, 1)
    [t, sol] = ode45(duffingEquations, time_span, initial_conditions(i, :));
    plot(sol(:, 1), sol(:, 2), 'Color', colors(i, :), 'LineStyle', line_styles{i}, 'LineWidth', 1);
end

% Customize the plot
axis([-2 2 -2 2]);
grid on;
grid minor;
set(gca, 'XTick', -2:0.5:2, 'YTick', -2:0.5:2, 'GridLineStyle', ':');
set(gca, 'Box', 'on'); 

% Add arrows
for i = 1:size(initial_conditions, 1)
    [t, sol] = ode45(duffingEquations, time_span, initial_conditions(i, :));
    % Adding arrows at specific points
    for j = 1:5:length(sol)-1
        % Scale factor for arrows
        arrow_scale = 0.1;
        % Compute the direction vectors
        u = sol(j+1, 1) - sol(j, 1);
        v = sol(j+1, 2) - sol(j, 2);
        % Normalize the direction vectors for consistent arrow size
        norm_factor = sqrt(u^2 + v^2);
        if norm_factor > 0
            u = arrow_scale * u / norm_factor;
            v = arrow_scale * v / norm_factor;
            quiver(sol(j, 1), sol(j, 2), u, v, 0, 'Color', colors(i, :), 'LineWidth', 1);
        end
    end
end

% Labels and title
xlabel('$U_n$', 'Interpreter', 'latex');
ylabel('$U_n^\prime$', 'Interpreter', 'latex'); 
title('Duffing Oscillator Phase Portrait');

hold off;

[Athanasios Paraskevopoulos]

@Sam Chak I would appreciate an example to understand annotation in a problem like this. I have read documentation and it is a bit different than my purpose

[Sam Chak]

Hi @Athanasios Paraskevopoulos

There are a few approaches to generating '→' in MATLAB, and it is common to use the 'quiver()' and the 'annotation()' functions, unless you wish to manually draw the arrows using the primitive 'line()' function.

The quiver() function is a straightforward way to create vector field-style arrows, where the coordinates of the arrowheads are specified using data points. While you can easily determine the start point of the arrow tail, you'll need calculate the direction of the arrow. It is also relatively difficult to adjust the scale of the arrowhead's size for short quivers.

Alternatively, the annotation() function can also do the job and it has a better-looking arrowhead, but as you may have noted, the coordinate system is normalized to the figure itself rather than the data points. This can make the placement a bit more challenging, as the coordinates range from (0,0) in the lower-left corner to (1,1) in the upper-right corner of the figure.

Here is an example where I placed a short arrow at the data point origin (x = 0, y = 0)

fig = figure;
panel = uipanel(fig);

x   = linspace(-4, 4, 8001);
y   = tanh(x);

axes('Parent', panel);

plot(x, y), grid on, grid minor

%% Arrow's coordinates that do not depend on Data Points
x_posA  =  0.1;
y_posA  = tanh(x_posA) - x_posA;
x_posB  =  x_posA + 0.1;
y_posB  = tanh(x_posB) - 0.1;

%% Calculate relative coordinates
x_relA  = (x_posA - min(x))/(max(x) - min(x));
y_relA  = (y_posA - min(y))/(max(y) - min(y));
x_relB  = (x_posB - min(x))/(max(x) - min(x));
y_relB  = (y_posB - min(y))/(max(y) - min(y));

% Add arrowhead annotation at x = 0, y = 0
annotation(panel, 'arrow', [x_relA, x_relB], [y_relA, y_relB], 'HeadStyle', 'vback2', 'HeadWidth', 10, 'HeadLength', 10, 'Color', 'red')

% axis tight
hold off

Answer 1

[Sam Chak]

Hi @Athanasios Paraskevopoulos

This method utilizes the 'arrow()' function, which was created by Professor @Erik Johnson. You can download his MATLAB m-file from the MATLAB Central File Exchange, and acknowledge his contribution.

It is also worth noting that it is unnecessary to solve the differential equation using ode45 solver, as the 'stream2()' function can be used to generate the 2D data points for plotting purposes. A demonstration for the Duffing oscillator is provided below. You may choose to optimize the code using a for-loop approach.

%% Parameters
omega_d_squared = 1;
beta    = 1;
delta   = 0.5;

%% Duffing equations
[X, Y]  = meshgrid(-2:4/14:2);
U       = Y; 
V       = omega_d_squared*X - beta*X.^3 - delta*Y;

hold on
lineobj = streamslice(X, Y, U, V);
set(lineobj, 'Color', '#A9E0FB');
% set(lineobj, 'Density', 0.5);

%% Trajectory #1
xStart  = -1.5;
yStart  = -1.5;

    %% generate 2D data of the vector field
    verts   = stream2(X, Y, U, V, xStart, yStart);
    array   = verts{1};
    x       = array(:,1);
    y       = array(:,2);
    
    %% plot from start position to the arrowtail location
    n       = 60;
    idx     = round(length(x)/n);
    idy     = round(length(y)/n);
    plot(x(1:idx), y(1:idy), 'linewidth', 1.5, 'Color', '#265EF5')
    
    %% plot arrow at the desired location
    warning('off')
    arrow([x(idx), y(idy)], [x(idx+5), y(idy+5)], 'linewidth', 1.5, 'Color', '#265EF5')
    
    %% plot from arrowhead location to the end location of the trajectory
    plot(x(idx+5:end), y(idy+5:end), 'linewidth', 1.5, 'Color', '#265EF5')

%% Trajectory #2
xStart  = -0.01;
yStart  = 0;

    %% generate 2D data of the vector field
    verts   = stream2(X, Y, U, V, xStart, yStart);
    array   = verts{1};
    x       = array(:,1);
    y       = array(:,2);
    
    %% plot from start position to the arrowtail location
    n       = 250;
    idx     = round(length(x)/n);
    idy     = round(length(y)/n);
    plot(x(1:idx), y(1:idy), 'linewidth', 1.5, 'Color', '#FA477A')
    
    %% plot arrow at the desired location
    warning('off')
    arrow([x(idx), y(idy)], [x(idx+5), y(idy+5)], 'linewidth', 1.5, 'Color', '#FA477A')
    
    %% plot from arrowhead location to the end location of the trajectory
    plot(x(idx+5:end), y(idy+5:end), 'linewidth', 1.5, 'Color', '#FA477A')

%% Trajectory #3
xStart  = 0.01;
yStart  = 0;

    %% generate 2D data of the vector field
    verts   = stream2(X, Y, U, V, xStart, yStart);
    array   = verts{1};
    x       = array(:,1);
    y       = array(:,2);
    
    %% plot from start position to the arrowtail location
    n       = 250;
    idx     = round(length(x)/n);
    idy     = round(length(y)/n);
    plot(x(1:idx), y(1:idy), 'linewidth', 1.5, 'Color', '#E6AC00')
    
    %% plot arrow at the desired location
    warning('off')
    arrow([x(idx), y(idy)], [x(idx+5), y(idy+5)], 'linewidth', 1.5, 'Color', '#E6AC00')
    
    %% plot from arrowhead location to the end location of the trajectory
    plot(x(idx+5:end), y(idy+5:end), 'linewidth', 1.5, 'Color', '#E6AC00')

%% Trajectory #4
xStart  = 1.5;
yStart  = 1.5;

    %% generate 2D data of the vector field
    verts   = stream2(X, Y, U, V, xStart, yStart);
    array   = verts{1};
    x       = array(:,1);
    y       = array(:,2);
    
    %% plot from start position to the arrowtail location
    n       = 60;
    idx     = round(length(x)/n);
    idy     = round(length(y)/n);
    plot(x(1:idx), y(1:idy), 'linewidth', 1.5, 'Color', '#23BB98')
    
    %% plot arrow at the desired location
    warning('off')
    arrow([x(idx), y(idy)], [x(idx+5), y(idy+5)], 'linewidth', 1.5, 'Color', '#23BB98')
    
    %% plot from arrowhead location to the end location of the trajectory
    plot(x(idx+5:end), y(idy+5:end), 'linewidth', 1.5, 'Color', '#23BB98')    

hold off
grid on, axis([-2 2 -2 2]), 

%% Labels
xlabel('$U_n$',        'Interpreter', 'latex')
ylabel('$U_n^\prime$', 'Interpreter', 'latex')
title('Phase Portrait of Duffing Oscillator')

[Athanasios Paraskevopoulos]

@Sam Chak Thank you for your explanation. It was very helpful