function two_variables_RK
clear all
close all
clc
format long
%比较Euler法、四阶RK算法-----

time = 10; N = 40000;
h = time/N;
t = 0:h:time;

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

%--------------解析解------------------------
m=1;    %谐振子的质量为1
k=4*pi^2; %弹簧的劲度系数
A=sqrt(y0^2+p0^2/(m*k)); %谐振子的振幅
pha=acos(y0/A);          %谐振子的初始相位，若x0与p0同号，pha为负，否则pha为正
y_e=A*cos(sqrt(k/m)*t+pha); %谐振子的精确解
p_e=-m*A*sqrt(k/m)*sin(sqrt(k/m)*t+pha); %动量的精确解
E_e=p_e.^2/(2*m)+0.5*k*y_e.^2; %能量的精确值
E_e(1)

%----------Euler算法------------------------------
yyy=zeros(1,N+1);
ppp=zeros(1,N+1);
yyy(1)=y0;
ppp(1)=p0;
for n=1:N %Euler法迭代公式
    yyy(n+1)=yyy(n)+h*ppp(n);
    ppp(n+1)=ppp(n)-h*4*pi^2*yyy(n);
end
E_euler=ppp.^2/(2*m)+0.5*k*yyy.^2; %能量的数值解
figure,subplot(2,2,1),plot(t,yyy,'r',t,y_e,'b')
title('谐振子的位置（Euler算法与精确解比较）')
legend('Euler算法','精确解')
subplot(2,2,2),plot(t,ppp,'r',t,p_e,'b')
title('谐振子的动量（Euler算法与精确解比较）')
legend('Euler算法','精确解')
subplot(2,2,3),plot(yyy,ppp,'r',y_e,p_e,'b')
title('谐振子的相图（Euler算法与精确解比较）')
legend('Euler算法','精确解')
subplot(2,2,4),plot(t,E_euler-E_e,'r')
title('误差')

%----------Euler算法的另一种程序写法，列向量形式------------------------------
f=@(t,y,p)[p; -4*pi^2*y]; % 定义 f，写成列向量形式
z0 = [y0; p0];      %初始条件，y0=5,p0=0
z_euler=zeros(2,N+1);
z_euler(:,1)=z0;
for n=1:N %Euler法迭代公式
    z_euler(:,n+1)=z_euler(:,n)+h*f(t(n),z_euler(1,n),z_euler(2,n));
end
y_euler = z_euler(1,:); 
p_euler = z_euler(2,:); 
E_euler_1=p_euler.^2/(2*m)+0.5*k*y_euler.^2; %能量的数值解
figure,subplot(2,2,1),plot(t,y_euler,'r',t,y_e,'b')
title('谐振子的位置（程序写法2）')
legend('Euler算法','精确解')
subplot(2,2,2),plot(t,p_euler,'r',t,p_e,'b')
title('谐振子的动量（程序写法2）')
legend('Euler算法','精确解')
subplot(2,2,3),plot(y_euler,p_euler,'r',y_e,p_e,'b')
title('谐振子的相图（程序写法2）')
legend('Euler算法','精确解')
subplot(2,2,4),plot(t,E_euler_1-E_e,'r')
title('误差')

figure;plot(t,yyy-y_euler,'r')


%--------------四阶RK算法------------------------
yy=zeros(1,N+1);
pp=zeros(1,N+1);
yy(1)=y0;
pp(1)=p0;
%f=@(t,y,p)[p; -4*pi^2*y]; % 定义 f，写成列向量形式
%z0 = [y0; p0];      
%初始条件，y0=5,p0=0，已经在Euler法的第二种写法中定义
z  = myrungekutta(f, t, z0);
yy = z(1,:); 
pp = z(2,:); 
E_4rk=pp.^2/(2*m)+0.5*k*yy.^2; %能量的数值解

figure,subplot(2,2,1),plot(t,yy,'r',t,y_e,'b')
title('谐振子的位置（四阶RK算法与精确解比较）')
legend('四阶RK算法','精确解')
subplot(2,2,2),plot(t,pp,'r',t,p_e,'b')
title('谐振子的动量（四阶RK算法与精确解比较）')
legend('四阶RK算法','精确解')
subplot(2,2,3),plot(yy,pp,'r',y_e,p_e,'b')
title('谐振子的相图（四阶RK算法与精确解比较）')
legend('四阶RK算法','精确解')
subplot(2,2,4),plot(t,E_4rk-E_e,'r')
title('误差')

figure,plot(t,yy-y_e,'b')
 
function z= myrungekutta(f, t, z0)
n = length(t);
h = t(2)-t(1);
z=zeros(2,n);
%定义一个2行n列的数组，2为未知数的个数，也即方程组中方程的个数
z(:,1)=z0;
for i = 1:n-1
    %f=@(t,y,p)[p; -4*pi^2*y]
    k1=f(t(i), z(1,i),z(2,i));
    k2=f(t(i)+0.5*h, z(1,i) + h*k1(1)/2, z(2,i) + h*k1(2)/2);
    k3=f(t(i)+0.5*h, z(1,i) + h*k2(1)/2, z(2,i) + h*k2(2)/2);
    k4=f(t(i)+h, z(1,i) + h*k3(1), z(2,i) + h*k3(2));
    z(:,i+1)=z(:,i) + h*(k1 + 2*k2 + 2*k3 + k4)/6;
    %z(1,i+1)=z(1,i) + h*(k1(1) + 2*k2(1) + 2*k3(1) + k4(1))/6;
    %z(2,i+1)=z(2,i) + h*(k1(2) + 2*k2(2) + 2*k3(2) + k4(2))/6;
end

















