Plotting Phase Portrait of Duffing Equation
23 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Athanasios Paraskevopoulos
el 22 de Mayo de 2024
Comentada: Athanasios Paraskevopoulos
el 25 de Mayo de 2024
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1701901/image.png)
1. Equation Details:
- The first equation:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1701906/image.png)
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1701911/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1701916/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1701921/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1701926/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1701931/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1701936/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1701941/image.png)
- Boundary conditions:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1701946/image.png)
These conditions specify that $ \Phi $ is zero at the boundaries of the domain.
2. Special Case (Non-coupled Particles,
:)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1701951/image.png)
- In the limit $ k = 0 $, the spatial derivative terms disappear, simplifying the equation to:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1701956/image.png)
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1701961/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1701966/image.png)
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1701976/image.png)
4 comentarios
Sam Chak
el 25 de Mayo de 2024
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
Respuesta aceptada
Sam Chak
el 25 de Mayo de 2024
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')
Más respuestas (0)
Ver también
Categorías
Más información sobre 2-D and 3-D Plots en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!