## Example of numerical solution

The Duffing equation is a model for mechanical oscillations with a nonlinear spring. The dynamics of this system, at a certain parameter is described by:

 (5.7)

To use MATLAB to compute the direction field, trajectories and solutions of this system, we rewrite this system in the state space form , where . The nonlinear function becomes

 (5.8)

This equation is implemented in the MATLAB routine duff.m, containing
function qdot=duffing(t,q);
qdot=[q(2);
-0.1*q(2)+q(1)-q(1)^3];

In the following script file, the above mentioned function file is used, to plot first the direction field, afterwards the trajectory over time, and finally the solution curves in phase plane, together with the direction field. Since is not explicitly dependent on , the script uses duffing(0,q) to compute the direction field, this is valid at all times.
%% script file to analyse duffing equation.

clear all;
close all;

%% create grid for direction plot
[Q1,Q2]=meshgrid(-1.5:.3:1.5,-0.7:.14:0.7);

%% create matrices with elements of vector Qd
Qd1=zeros(size(Q1));
Qd2=zeros(size(Q2));
for i=1:11;
for j=1:11;
Qdot=duffing(0,[Q1(i,j);Q2(i,j)]);
Qd1(i,j)=Qdot(1);
Qd2(i,j)=Qdot(2);
end
end

%% plot direction plot
figure(1);
quiver(Q1,Q2,Qd1,Qd2);
xlabel('q_1');
ylabel('q_2')
title('Direction field of Duffing equation');

%% create time vector and numerical solution with ode45
[t, Q]=ode45('duffing',[0,20],[1.5;0]);

%% plot trajectory in time
figure(2)
plot(t,Q(:,1))
xlabel('t');
ylabel('q_1');

figure(3)
plot(t,Q(:,2))
xlabel('t');
ylabel('q_2');

%% plot trajectory in phase space
figure(4);
quiver(Q1,Q2,Qd1,Qd2);
xlabel('q_1=x');
ylabel('q_2=dotx')
title('Solution curve in phase plane');
hold on;
plot(Q(:,1),Q(:,2),'k');


Running these files yields the Figures 5.4-5.6.

Esteur 2010-03-22