function jyjs_2
clear 
clc
close all
format long

options=odeset('RelTol',1e-14);
[X,Y]=ode45(@yjs,[0,10],[5;0],options);

%--------------解析解------------------------
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为正
t=X';
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_m=Y(:,2)'.^2/(2*m)+0.5*k*Y(:,1)'.^2; %能量的数值解
figure;
subplot(2,2,1);plot(X,Y(:,1),'r',X,y_e,'b')
subplot(2,2,2);plot(X,Y(:,2),'r',X,p_e,'b')
subplot(2,2,3);plot(Y(:,1),Y(:,2))
subplot(2,2,4);plot(X',E_m-E_e,'r')

%定义句柄函数
function ydot=yjs(t,y)
ydot=[y(2);-4*pi^2*y(1)];
