Forced damped Pendulum

First, for a samplling of values of the parameter A, we plot simulations. We take samples of the solution once per period of the forcing using the events location feature of ode45.

b=1/2;

omega=2/3;

A_vector=[0.9 1.07 1.15 1.35 1.46 1.5];

behaviors=["Periodic orbit","Two-cycle","Chaos","Periodic Orbit","Four-cycle","Chaos"];

nA=length(A_vector);

y0=[0 1];

tSpan=[0 1000];

tTransient=200;

for k=1:nA

A=A_vector(k);

f=@(t,y)pendForceDamp(t,y,b,A,omega);

fEvent=@(t,y)pendulumEvent(t,y,omega);

options=odeset('Events',fEvent,'RelTol',1e-8,'AbsTol',1e-8);

[t,y,te,ye]=ode45(f,tSpan,y0,options);

% Throw away a transient

spot=find(t>tTransient,1,'first');

y=y(spot:end,:);

% Throw away a transient from the events time series

spot=find(te>tTransient,1,'first');

ye=ye(spot:end,:);

figure

plot(y(:,1),y(:,2),ye(:,1),ye(:,2),'o')

title(sprintf('$A=%0.2f$: %s',A),behaviors(k))

xlabel('$\theta$');ylabel('$\frac{d\theta}{dt}$')

figure

plot(mod(ye(:,1),2*pi),ye(:,2),'o')

axis([0 2*pi min(y(:,2)) max(y(:,2))])

title(sprintf('$A=%0.2f$: %s',A),behaviors(k))

xlabel('$\theta$');ylabel('$\frac{d\theta}{dt}$')

y0=y(end,:); % use the final solution as the new initial condtion

end

Bifurcation Diagram

Here we see some period-doubling but we do not see a diagram equivalent to the logistic map.

nA=301;

a0=0.9;

a1=1.5;

A_vector=linspace(a0,a1,nA);

tSpan=[0 2000];

tTransient=1000;

figure;clf;hold on;

for k=1:nA

A=A_vector(k);

f=@(t,y)pendForceDamp(t,y,b,A,omega);

fEvent=@(t,y)pendulumEvent(t,y,omega);

options=odeset('Events',fEvent,'RelTol',1e-8,'AbsTol',1e-8);

[~,~,te,ye]=ode45(f,tSpan,y0,options);

% Throw away a transient from the events time series

spot=find(te>tTransient,1,'first');

ye=ye(spot:end,:);

plot(ones(length(ye),1)*A,ye(:,2),'b.')

end

xlabel('$A$');ylabel('$\frac{d\theta}{dt}$')

Functions called

function dydt=pendForceDamp(t,y,b,A,omega)
dydt=zeros(2,1);
dydt(1) = y(2);
dydt(2) = -b* y(2) - sin(y(1)) + A*sin(omega*t);
end
function  [value, isterminal,direction]=pendulumEvent(t,~,omega)
value= sin(omega*t);
isterminal = 0;
direction = 1;
end