清空筛选条件
清空筛选条件

Solving Duffing Equation with the new framework for ODEs

20 次查看（过去 30 天）
Athanasios Paraskevopoulos
Athanasios Paraskevopoulos 2024-8-17，14:34
评论： Athanasios Paraskevopoulos 2024-8-17，16:09
Ran in:
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;

请先登录，再回答此问题。

回答（1 个）

Steven Lord
Steven Lord 2024-8-17，15:41
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?
  1 个评论
Athanasios Paraskevopoulos
Athanasios Paraskevopoulos 2024-8-17，16:09
Ran in:
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by