# Plotting Phase Portrait of Duffing Equation

75 views (last 30 days)
Athanasios Paraskevopoulos on 22 May 2024 at 22:55
Commented: Athanasios Paraskevopoulos on 25 May 2024 at 5:05
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
##### 4 CommentsShow 2 older commentsHide 2 older comments
Athanasios Paraskevopoulos on 24 May 2024 at 23:34
@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 on 25 May 2024 at 1:30
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

Sign in to comment.

### Accepted Answer

Sam Chak on 25 May 2024 at 3:46
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')
##### 1 CommentShow -1 older commentsHide -1 older comments
Athanasios Paraskevopoulos on 25 May 2024 at 5:05
@Sam Chak Thank you for your explanation. It was very helpful

Sign in to comment.

### Categories

Find more on 2-D and 3-D Plots in Help Center and File Exchange

R2024a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!