test1.m

R=1; % pendulum length
g=9.8; % gravity
al=0.5; % friction

% initial condition:
theta=2*pi/3;
dtheta=0;

% base:
fb=113; % frequency
ab=0.05; % amplitude
x=@(t) 0; % x
y=@(t) ab*sin(2*pi*fb*t); % y
ddx=@(t) 0; % x- acceleration
ddy=@(t) -ab*(2*pi*fb)^2*sin(2*pi*fb*t); % y-acceleration

dt=0.00013; % time step
tm=8; % end time

% mass position:
xm=x(0)+R*cos(theta);
ym=y(0)+R*sin(theta);

% inital plot
hp=plot([x(0) xm],[y(0) ym],'or-');
axis equal;
xlim(1.3*[-R R]);
ylim(1.3*[-R R]);

% add arrow:
hold on;
rr=[R R];
rr1=rr+[0,-0.1*R];
rr2=rr+[0,-0.1*R]+[-0.03*R 0.03*R];
rr3=rr+[0,-0.1*R]+[0.03*R 0.03*R];
plot([rr(1) rr1(1)  rr2(1)],[rr(2) rr1(2) rr2(2)],'k-');
plot([rr1(1)  rr3(1)],[rr1(2) rr3(2)],'k-');
rr3=(rr1+rr2)/2+[0.05*R 0];
 % g-symbol:
text(rr3(1),rr3(2),'g','horizontalalignment','left');


lc=1; % loop counter
dlc=10; % how many frames to skeep
for t=0:dt:tm
    rs=(-g*sin(theta)-ddx(t)*cos(theta)-ddy(t)*sin(theta))/R-al*dtheta; % right side
    
    % integrate (Euler method):
    dtheta=dtheta+dt*rs;
    theta=theta+dt*dtheta;
    
    % mass position:
    xm=x(t)+R*sin(theta);
    ym=y(t)-R*cos(theta);
    
    % display:
    if mod(lc,dlc)==0 % skeep frames
        set(hp,'XData',[x(t) xm],'YData',[y(t) ym]);
        xlim(1.3*[-R R]);
        ylim(1.3*[-R R]);
        drawnow;
    end
    
    lc=lc+1;
end