Example of numerical solution

(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

(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

`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