# Nonlinear ODE to solve Duffing's equation

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

• Your code looks more or less correct but you have omitted so many details it is difficult to say more about where the problem might be. I don't see how you are passing the values of the parameters (e.g. k, c) into ode_fun, for example. I think the actual values of those parameters (which you don't supply) may also give some idea about why ode45 thinks the solution is linear over the range t=0,T. – Bill Greene Nov 15 '14 at 16:57
• I'll give you also the numbers I used. – Rhei Nov 15 '14 at 17:40
• I think those numbers are OK. What version of MATLAB are you using? I think you will need to post your complete code in runnable form for someone to be able to help you further. – Bill Greene Nov 15 '14 at 18:57
• anyone who knows the code for frequency responce of the Duffing equation to get the plot like the one on scholar pedia – M91 Sep 18 '18 at 13:10

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.