ex22ch2.m

these codes are for solving OED with matlab

function ex22ch2
global f2 Omega2 Omega3 zeta f3
f2 = 0.05;
Omega2 = 1;
Omega3 = 2;
zeta = 0.1;
f3s = [0.50; 0.92; 0.99; 1.005; 1.35];
xp0s = [-0.60; -0.80; -0.59; -0.80; 0];
options = odeset('RelTol',1e-8,'AbsTol',1e-8);
for i = 1:5
   f3 = f3s(i); 
   [t,x] = ode45(@odes,[0, 500],[0; xp0s(i)],options);
   ndx = find(t >= 400);
   figure
   plot(x(ndx,1),x(ndx,2))
   xlabel('x(t)')
   ylabel('x''(t)')
   title(['System with two frequency excitation, ',...
          'f_3 = ',num2str(f3),'.'])
   drawnow
end
%============================================================
function dydt = odes(t,y)
global f2 Omega2 Omega3 zeta f3
dydt = [y(2); 0];
dydt(2) = -2*zeta*y(2) - y(1) - 1.5*y(1)^2 - 0.5*y(1)^3 + ...
           f2*cos(Omega2*t) + f3*cos(Omega3*t);