function dirven_pendulum
clear all
close all
clc

q=0; 
b=0; 
omega0=2/3; %周期性强迫外力的频率
theta0=pi*3/180; %初始角度3度，以弧度为单位
g=9.8;
l=1;

time = 20; N = 20000;
h = time/N;
t = 0:h:time;

f=@(t,x)[x(2), -q*x(2)-g/l*sin(x(1))+b*cos(omega0*t)];%定义句柄函数
z0 = [theta0, 0]; %初始条件，注意角度采用弧度单位，角速度为0

z = myrungekutta(f, t, z0);
theta = z(:,1);
omega = z(:,2); 

subplot(2,1,1) 
plot(t, theta,'b',t,omega,'r') 
title('驱动单摆的角度与角速度随时间的变化关系')
legend('角度','角速度')
subplot(2,1,2)
plot(theta, omega)
%plot(theta(N/2:N), omega(N/2: N))
title('角度-角速度相图')
figure;
subplot(2,1,1),plot(t,theta','b',t,theta0*cos(sqrt(g/l)*t),'r')
title('角度演化的数值与解析比较')
legend('数值结果','解析结果')
subplot(2,1,2),plot(t,theta'-theta0*cos(sqrt(g/l)*t),'r')

function z = myrungekutta(f, t, z0)
n = length(t);
h = t(2)-t(1);
z(1,:)=z0;
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

  
