function dirven_pendulum
clear all
close all
clc

q=0; 
b=0; 
w0=2/3; %周期性强迫外力的频率
theta0=linspace(pi/360,pi-pi/360,30); %初始角度，以弧度为单位
omega0=0; %自由释放，角速度为0

time = 20; N = 50000;
h = time/N;
t = 0:h:time;
f=@(t,x)[x(2), -q*x(2)-sin(x(1))+b*cos(w0*t)];
%定义句柄函数

for i=1:length(theta0)
    z0 = [theta0(i), omega0]; %初始条件，注意角度采用弧度单位，角速度为0
    z  = myrungekutta(f, t, z0);
    theta = z(:,1)';
    omega = z(:,2)';
    for j=1:N-1
        if omega(j+1)*omega(j)<0 %符号相反表示钟摆开始改变运动方向
            %该程序的缺点是RK算法会计算所有时间，不论钟摆完成多少个周期运动
           TT(i)=(2*j+1)*h
           break;  %跳出当前所在for循环
        end
    end
end
plot(theta0,TT,'o-');
title('驱动单摆的初始角度与周期的关系')
xlabel('初始角度'); ylabel('周期');


function z = myrungekutta(f, t, z)
n = length(t);
h = t(2)-t(1);
for i = 1:n-1
    k1=f(t(i), z(i,:));
    k2=f(t(i)+0.5*h, z(i,:) + h*k1/2);
    k3=f(t(i)+0.5*h, z(i,:) + h*k2/2);
    k4=f(t(i)+h, z(i,:) + h*k3);
    z(i+1,:)=z(i,:) + h*(k1 + 2*k2 + 2*k3 + k4)/6;
end

  
