%利用四阶RK算法求解含阻尼和周期性驱动力下的
%周期为1的一维谐振子运动
function two_variables_RK
clear all
close all
clc

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

y0=10;   %初始位置,可以任意设置
p0=20;   %初始动量,可以任意设置
yy=zeros(1,N+1);
pp=zeros(1,N+1);

b=500;   %周期性驱动力的振幅
w=pi/10; %周期性驱动力的频率，振动周期为10
f=@(t,x)[x(2), -4*pi^2*x(1)-0.6*x(2)+b*cos(w*t)]; 
% 定义 f,其中-0.6*x(2)为阻尼项，b*cos(w*t)为驱动力项
z0 = [y0, p0]; 
z = myrungekutta(f, t, z0);
yy = z(:,1)'; 
pp = z(:,2)'; 

subplot(2,1,1),plot(t,yy,'r')
title('y(t) line')
subplot(2,1,2),plot(t,pp,'r')
title('p(t) line')
figure,plot(yy,pp,'r')
figure,comet(yy,pp)

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