Nonlinear ODE to solve Duffing's equation

Question

I am trying to solve the Duffing's equation in MATLAB.

$ m\ddot{y}+c\dot{y}+ky+k_{3}y^{3} = f(t) $

where $ f(t) = A \sin{\omega t}$

To do that I wrote a function to be given to the ode45.

function [F] = ode_fun(t,y)

F = zeros(size(y));

F(1) = y(2);
F(2) = -k*y(1)-c*y(2)-k3*y(1)^3+A*sin(omega*t);

end

And then I call the ode function:

[t,u] = ode45('ode_fun',[0 T],[0 0]);

Is this the right way to do it?

Because I am expecting to obtain a cubic function when plotting the "force vs displacement" plot but I actually obtain a straight line.

EDIT: I'll give you also the numerical values I used for the simulation.

Basically I started from a linear equation I found in an example and then I added the nonlinear term to check if my code was able to distinguish between linear and nonlinear dynamic systems.

The values of the linear case are: $m = 1$ , $k = 10^4$, $ c = 40$.

Then I added the nonlinear term as $ k_3 = 11*10^4$

As forcing term I decided to use a sine wave having $A = 30$ and $\omega = 10$.

The final simulation time is $T = 60$

EDIT: the complete code I am testing is the one I talked about in this other question I posted 2D cross section from 3D surface

I am trying to understand if this procedure works for nonlinear parameter identification

Answer 1

Cubics can look quite linear in certain regimes. Your linear stiffness parameters, $k$ and $k_3$, are rather high. If you look at the the range of the position state variable with your specified values, you'll notice that the system is only deviating from zero by about $\pm 3\times10^{-3}$. Physically, this is effectively in the linear regime of the system so any nonlinear components will be very small. Reducing $k$ by an order of magnitude, the oscillations are larger and the plot of force (acceleration) vs. position starts to look cubic after the system settles into steady-state:

Some suggestions.
1) I may be a good idea to start using a stiff solver for a this system, depending on what your parameters are. In Matlab, ode15s is a good choice.

2) Looking at your code in your other question, you really should be passing your parameter values as variable rather than hard coding them in multiple places. See this article on parametrizing functions in Matlab for details. For your code:

tspan = [0 60];
k = 1e3;
c = 40;
k3 = 1.1e5;
A = 30;
omega = 10;

ode_fun = @(t,u,k,c,k3,A,omega)[u(2);-k*u(1)-c*u(2)-k3*u(1).^3+A*sin(omega*t)];

[t,u] = ode15s(@(t,u)ode_fun(t,u,k,c,k3,A,omega),tspan,[0 0]);

figure;
plot(t,u);
xlabel('Time');
ylabel('State');

acc = -k*u(:,1)-c*u(:,2)-k3*u(:,1).^3+A*sin(omega*t);
figure;
plot(u(:,1),acc);
xlabel('Position');
ylabel('Acceleration');

Here, I've used the the ODE function to directly obtain the acceleration.

3) You may want to specify your integration times (tspan) as a vector with more than two elements that you can ensure a particular output step size. (N.B., this not the actual step size used internally by the integrator – see here.) This will help if you use a finite difference scheme for calculating accelerations and it can produce smoother plots in some cases.