How plot the systems with perturbation term?

is been a while i am searching for find a good option of ploting in my search i found some usefull one but i need to seperate them becuase background size of picture is not what i want and i want them solo and if there is any other better option of plotting the system please provide it here thanks

de1 = diff(u(t), t) == v(t);
de2 = diff(v(t), t) == -K(1) * u(t) ^ 3 + K(2) * u(t) + epsilon * sin(theta * t);

[Torsten]

I don't know the meaning of the zeta and gamma parameters.

[salim]

@Torsten

% Chaotic System Visualization - Forced Duffing Oscillator
% Based on: du/dt = v, dv/dt = -K1*u^3 + K2*u + epsilon*sin(theta*t)
%% Parameters (adjust these to explore different chaotic behaviors)
K1 = 1.0;      % Cubic nonlinearity coefficient
K2 = -1.0;     % Linear restoring force coefficient  
epsilon = 1.0;  % Forcing amplitude
theta = 1.2;   % Forcing frequency

[Torsten]

zeta_1, zeta_2, zeta_3, gamma_1, gamma_2, gamma_3 - what do they mean ?

[salim]

[Torsten]

Something like this ?

% Chaotic System Visualization - Forced Duffing Oscillator
% Based on: du/dt = v, dv/dt = -K1*u^3 + K2*u + epsilon*sin(theta*t)
%% Parameters (adjust these to explore different chaotic behaviors)
K1 = 1.0;      % Cubic nonlinearity coefficient
K2 = -1.0;     % Linear restoring force coefficient  
epsilon = 1.0;  % Forcing amplitude
theta = 1.2;   % Forcing frequency
f = @(t,y)[y(2);-K1*y(1)^3+K2*y(1)+epsilon*sin(theta*t)];
tspan = [0 100];
y0 = [1;0];
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
[T,Y] = ode45(f,tspan,y0);
plot(Y(:,1),Y(:,2))

●

Answer 1

[Sulaymon Eshkabilov]

You can work around what Torsten suggested to enhance the resolution of simulation results and get the simulation for other values of epsilon, e.g.:

% Chaotic System Visualization - Forced Duffing Oscillator
% Based on: du/dt = v, dv/dt = -K1*u^3 + K2*u + epsilon*sin(theta*t)
%% Parameters (adjust these to explore different chaotic behaviors)
K1 = 1.0;      % Cubic nonlinearity coefficient
K2 = -1.0;     % Linear restoring force coefficient  
epsilon = 0.5:.5:2;  % Forcing amplitude
theta = 1.2;   % Forcing frequency
LT = {'r-'; 'g-'; 'b-'; 'm-'}; % Line color spec
for ii  = 1:numel(epsilon)
f = @(t,y)[y(2);-K1*y(1)^3+K2*y(1)+epsilon(ii)*sin(theta*t)];
tspan = linspace(0,100, 1e5);
y0 = [1;0];
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
[T,Y] = ode45(f,tspan,y0);
figure(ii)
plot(Y(:,1),Y(:,2), LT{ii})
title(['Phase portrait: (\epsilon =  ' num2str(epsilon(ii)) ')'])
grid on
xlabel('u')
ylabel('v')
end

●

[salim]

@Sulaymon Eshkabilov thank you so much, can we have 3D chaotic plot of each the plot too ?

also do you have any idea about poincare map of the system? and if have a good shape of phase portrait it will be great i have it but i need more option of plot

Answer 2

[Sulaymon Eshkabilov]

Here is how 3D plot can be generated, e.g.:

% Chaotic System Visualization - Forced Duffing Oscillator
% Based on: du/dt = v, dv/dt = -K1*u^3 + K2*u + epsilon*sin(theta*t)
%% Parameters (adjust these to explore different chaotic behaviors)
K1 = 1.0;      % Cubic nonlinearity coefficient
K2 = -1.0;     % Linear restoring force coefficient  
epsilon = 0.5:.5:2;  % Forcing amplitude
theta = 1.2;   % Forcing frequency
LT = {'r-'; 'g-'; 'b-'; 'm-'}; % Line color spec
for ii  = 1:numel(epsilon)
f = @(t,y)[y(2);-K1*y(1)^3+K2*y(1)+epsilon(ii)*sin(theta*t)];
tspan = linspace(0,100, 1e5);
y0 = [1;0];
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
[Time,Y] = ode45(f,tspan,y0);
figure(ii)
plot3(Time, Y(:,1),Y(:,2), LT{ii})
title(['Phase portrait: (\epsilon =  ' num2str(epsilon(ii)) ')'])
grid on
xlabel('Time')
ylabel('u')
zlabel('v')
end

●

[Sulaymon Eshkabilov]

Here are some starting points on a Poincare map:

K1 = 1.0;      % Cubic nonlinearity coefficient
K2 = -1.0;     % Linear restoring force coefficient  
epsilon = 0.5:.5:2;  % Forcing amplitude
theta = 1.2;   % Forcing frequency
LT = {'ro'; 'g*'; 'bd'; 'mp'}; % Marker color for Poincare maps

T = 2*pi/theta;        % Forcing period
nPeriods = 300;        % Number of stroboscopic samples (adjust for clarity)

for ii  = 1:numel(epsilon)
    f = @(t,y)[y(2); -K1*y(1)^3 + K2*y(1) + epsilon(ii)*sin(theta*t)];
    tspan = linspace(0, nPeriods*T, 1e5);  % Simulate multiple periods
    y0 = [1; 0];
    options = odeset('RelTol',1e-10,'AbsTol',1e-10);
    [Time,Y] = ode45(f, tspan, y0, options);
    
    % Poincare section at multiples of forcing period T:
    poincare_points = [];
    t_samples = (0:nPeriods-1) * T;
    
    for k = 1:length(t_samples)
        % Find closest time index:
        [~, idx] = min(abs(Time - t_samples(k)));
        poincare_points = [poincare_points; Y(idx,1), Y(idx,2)];
    end
    
    % Plotting Poincare map simulation results:
    figure(ii)
    plot(poincare_points(:,1), poincare_points(:,2), LT{ii}, 'MarkerSize', 6)
    title(['Poincaré Map (\epsilon = ' num2str(epsilon(ii)) ')'])
    xlabel('u'); ylabel('v');
    axis tight
    grid on
end

●

[salim]

@Sulaymon Eshkabilov why our poincare is so different i find it from another paper the same problem but poincare is so different, what is make my poincare be like that ?