Solving Duffing Equation with the new framework for ODEs

I solved the duffing equation using the new framework for ODEs. Below is the code.

It works fine and plots as expected. However, I was wondering if there is a way to create a phase portrait in the interval t=[0 3000] without looking like a smudge.

% Define the parameters
delta = 0.1;
alpha = -1;
beta = 1;
gamma = 0.35;
omega = 1.4;

% Define the ODE function for the Duffing equation
duffingODE = @(t, y) [y(2); 
                      -delta*y(2) - alpha*y(1) - beta*y(1)^3 + gamma*cos(omega*t)];

% Initial conditions: [x(0), dx/dt(0)]
initialConditions = [0; 0];

% Create an ode object
F = ode(ODEFcn=duffingODE, InitialTime=0, InitialValue=initialConditions);

% Solve the equation over the interval [0, 3000]
sol = solve(F, 0, 3000);

% Interpolate the solution to get more points for plotting
timeFine = linspace(0, 300, 10000); % Create a fine time vector with 10,000 points
solutionFine = interp1(sol.Time, sol.Solution', timeFine)'; % Interpolate solution

% Plot the interpolated time series solution
figure;
subplot(2, 1, 1);
plot(timeFine, solutionFine(1, :), 'LineWidth', 1.5);
xlabel('Time');
ylabel('Displacement');
title('Interpolated Solution of the Duffing Equation');
grid on;

% Plot the interpolated phase portrait
subplot(2, 1, 2);
plot(solutionFine(1, :), solutionFine(2, :), 'LineWidth', 1.5);
xlabel('Displacement x(t)');
ylabel('Velocity dx/dt');
title('Interpolated Phase Portrait of the Duffing Equation');
grid on;

Answer 1

[Steven Lord]

FYI rather than interpolating the solution yourself with this code:

% Interpolate the solution to get more points for plotting
timeFine = linspace(0, 300, 10000); % Create a fine time vector with 10,000 points
solutionFine = interp1(sol.Time, sol.Solution', timeFine)'; % Interpolate solution

you can call solve on the ODE object with the timeFine vector directly, rather than just the start and end times. From the solve documentation page: "S = solve(F,t) computes the solution for the ODE represented by F at the specified time values in the vector t."

But to your main question: However, I was wondering if there is a way to create a phase portrait in the interval t=[0 3000] without looking like a smudge.

that's just a plain old line. Do you have a picture that shows the type of plot that you want to see, one that doesn't look "like a smudge"? Do one of the thumbnails on this documentation page look more like what you want?

[Athanasios Paraskevopoulos]

Hello @Steven Lord. First of all, thank you for your response.

When I don't use

% Interpolate the solution to get more points for plotting
timeFine = linspace(0, 300, 10000); % Create a fine time vector with 10,000 points
solutionFine = interp1(sol.Time, sol.Solution', timeFine)'; % Interpolate solution

The result is the following. Is there a way to have a smoother and clearer phase-plane? Or do I have to choose a specific interval ?

% Define the parameters
delta = 0.1;
alpha = -1;
beta = 1;
gamma = 0.35;
omega = 1.4;

% Define the ODE function for the Duffing equation
duffingODE = @(t, y) [y(2); 
                      -delta*y(2) - alpha*y(1) - beta*y(1)^3 + gamma*cos(omega*t)];

% Initial conditions: [x(0), dx/dt(0)]
initialConditions = [0; 0];

% Create an ode object
F = ode(ODEFcn=duffingODE, InitialTime=0, InitialValue=initialConditions);

% Solve the equation over the interval [0, 3000]
sol = solve(F, 0, 3000);

% Plot the time series solution
figure;
subplot(2, 1, 1);
plot(sol.Time, sol.Solution(1, :), 'LineWidth', 1.5);
xlabel('Time');
ylabel('Displacement');
title('Solution of the Duffing Equation');
grid on; % Add grid for better readability

% Plot the phase portrait
subplot(2, 1, 2);
plot(sol.Solution(1, :), sol.Solution(2, :), 'LineWidth', 1.5);
xlabel('Displacement x(t)');
ylabel('Velocity dx/dt');
title('Phase Portrait of the Duffing Equation');
grid on; % Add grid for better readability