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

time = 50; N = 10000;
h = time/N;
t = 0:h:time;

y0=2;   %初始位置,可以任意设置
p0=2;   %初始动量,可以任意设置

yy=zeros(1,N+1);
pp=zeros(1,N+1);

b=50;   %周期性驱动力的振幅
w=pi/5; %周期性驱动力的频率

f=@(t,x)[x(2), -4*pi^2*x(1)+b*cos(w*t)]; % 定义 f,其中b*cos(w*t)为驱动力项
z0 = [y0, p0]; %初始条件

z = myrungekutta(f, t, z0);
yy = z(:,1); 
pp = z(:,2); 

figure,plot(t,yy,'r')
title('坐标随时间的变化关系')
axis([0,50,-5,5])
figure,plot(t,pp,'r')
title('动量随时间的变化关系')
axis([0,50,-9,9])
figure,plot(yy,pp,'r')
title('坐标-动量相图')

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
