function sol = prob_4_8 % % Written by: % -- % John L. Weatherwax 2005-05-07 % % email: wax@alum.mit.edu % % Please send comments and especially bug reports to the % above email address. % %----- state = +1; options = ddeset( 'RelTol', 1e-5, 'AbsTol', 1e-8, 'Jumps', [ -10^(-6) ], 'Events', @events ); h6 = 10; sol = dde23(@ddes,[0.5 1], @history , [0, 60], options, state, h6 ); while( sol.x < 60 ) if( sol.ie(end) == 1 ) % Integration ended because of an event state = -state; sol = dde23(@ddes,[0.5 1], sol, [sol.x(end), 60], options, state, h6 ); end end yplot = sol.y; yplot(1,:) = (1e4)*yplot(1,:); yplot(2,:) = 0.5*yplot(2,:); yplot(4,:) = 10*yplot(4,:); figure; plot(sol.x,yplot,'-o'); grid on; title( [ 'h6 = ',num2str(h6) ] ); h6=300; sol = dde23(@ddes,[0.5 1], @history , [0, 60], options, state, h6 ); while( sol.x < 60 ) if( sol.ie(end) == 1 ) % Integration ended because of an event state = -state; sol = dde23(@ddes,[0.5 1], sol, [sol.x(end), 60], options, state, h6 ); end end yplot = sol.y; yplot(1,:) = (1e4)*yplot(1,:); yplot(2,:) = 0.5*yplot(2,:); yplot(4,:) = 10*yplot(4,:); figure; plot(sol.x,yplot,'-o'); grid on; title( [ 'h6 = ',num2str(h6) ] ); function v = history(t, state, h6 ) v = [ max(0,t+10^(-6)); 1; 1; 0 ]; function [value,isterminal,direction] = events(t,y,Z,state,h6) value = [ y(4)-0.1 ]; isterminal = 1; direction = 0; function dydt = ddes(t,y,Z,state,h6) h1 = 2; h2 = 0.8; h3 = 10^4; h4 = 0.17; h5 = 0.5; h7 = 0.12; h8 = 8; if( state == +1 ) xhi = 1; else xhi = (10/9)*(1-y(4)); end dydt = [ (h1-h2*y(3))*y(1); xhi*h3*Z(3,1)*Z(1,1) - h5*(y(2)-1); h4*(y(2)-y(3)) - h8*y(3)*y(1); h6*y(1) - h7*y(4); ];